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ABSTRACT 

World  leaders  depend  increasingly  on  remote  sensing  technologies  to  gather  information 
required  for  objective  decision  making  as  technologies  and  historical  events  blur  national 
boundaries.  The  Open  Skies  Treaty  depends  on  remote  sensing  in  the  form  of  aerial  photography 
to  enforce  its  bylaws  and  further  inspire  trust  and  cooperation  among  the  international  signatories. 
The  Open  Skies  Treaty  provides  guidelines  allowing  participants  to  fly  in  air  space  over  other 
participants’  countries  to  monitor  strategic  military  placement  and  development.  The  treaty 
restricts  the  ground  size  of  the  smallest  detail  recorded  by  these  aerial  imaging  systems  to  any  size 
larger  than  30  cm.  This  restriction  is  enforced  by  placing  a  lower  limit  on  the  altitude  at  which  a 
participating  aircraft  can  fly  and  it  is  computed  as  the  value  of  Current  techniques  rely  on 
human  photographic  interpreters  to  select  the  value  of  Hmm  for  every  calibration  pass  and  is  very 
resource  intensive.  The  Open  Skies  participants  are  investigating  machine  based  techniques  to 
supplement  the  traditional  human  role  in  an  effort  to  increase  the  objectiveness  of  the 
measurement. 

This  thesis  presents  a  software  tool  called,  ADiM,  a  man-in-the-loop,  algorithm  which 
manipulates  image  statistics  to  identify  the  orientation  and  width  of  individual  target  bar  groups 
from  digitized  images  of  aerial  photographs  of  Open  Skies  Treaty  calibration  triple  bar  target. 
ADiM  Hmin  results  achieved  an  88.6  percent  correlation  with  the  Open  Skies  Media  Processing 
Facility’s  computations. 
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Automatic  Digital  Processing  for  Calibration  Data  of  Open  Skies 

Treaty  Sensors 


I.  Introduction 

NASA's  Voyager  Space  Probes  I  and  II  possessed  the  ability  to  record  visual  imagery 
data  of  planets  Jupiter,  Uranus  and  their  moons  and  to  transmit  this  information  back  to  Earth. 
The  technologies  involved  in  unmanned  space  exploration,  such  as  the  Voyagers,  relied 
heavily  on  the  application  of  remote  sensing’  systems  [Wilhelms,  1984:  167;  Merchant, 

1996],  In  earth  orbit,  satellites  record  information  and  transmit  it  to  ground  based  receiving 
stations  using  the  same  remote  sensing  technology  as  Voyager.  Governments  and  industry 
both  benefit  from  these  remote  sensing  satellites  which  are  used  to  map  weather  systems  for 
predicting  weather,  or  to  map  large  areas  of  land  on  Earth  for  land  development  use  or 
environmental  measurements  [Merchant,  1996].  Similar  technology  can  be  used  to  collect 
imagery  of  Earth-based  regions  of  interest  by  the  United  States  and  many  other  countries  to 
assess  the  capability  and  intentions  of  potential  adversaries.  Remote  sensing  systems  are 
placed  on  aircraft  as  well  as  in  earth  orbit  and  can  )deld  similar  types  of  information. 

Aircraft  or  spacecraft  mounted  remote  sensing  systems  also  play  important  roles  in 
international  treaty  verification  [Merchant,  1996].  The  usefulness  of  aerial  imagery  in  the 

*  In  the  most  general  tenns  remote  sensing  is  “the  acquisition  of  information  about  an  object  without  physical 
contact”  [Simonett,  1983;  1]. 
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international  arena  is  accentuated  by  the  United  States’  National  Security  Policy  which  defines 
nine  fundamental  roles  for  the  United  States  Air  Force  (USAF).  As  one  of  these  nine 
fundamental  missions.  Air  Surveillance  and  Recoimaissance  requires  the  collection  and 
processing  of  information  from  airborne,  orbital,  and  surface  based  systems  to  provide 
strategic  and  tactical  intelligence  during  peace  and  war  [AFM  1-1,  1984;  3-5],  Since  1969, 
following  the  initial  Cold  War  Strategic  Arms  Limitations  Talks  (SALT  I)  and  the  Anti- 
Ballistic  Missile  Treaty  (ABM  Treaty)  between  the  US  and  the  former  Union  of  Soviet 
Socialist  Republics  (USSR),  remote  sensing  has  been  an  internationally  accepted  method  of 
treaty  verification  [Luttwak,  1996], 

1.1.  Open  Skies  Treaty 

In  1992,  the  United  States  of  America  became  a  signatory  to  the  Open  Skies  Treaty 
(OST).  OST  is  an  international  treaty  whose  provisions  grant  participating  countries^  the 
authority  to  collect  an  annual  quota  of  aerial  imagery  over  strategic  areas  within  the 
boundaries  of  other  signatory  countries.  The  OST  participants  aspire  to  improve  confidence 
between  their  respective  countries  [AF  News  Service,  1995],  The  following  text  is  an  excerpt 
from  the  treaty  and  reflects  its  ideals  and  purpose. 

Treaty  on  Open  Skies: 

The  States  concluding  this  Treaty 


2  This  list  of  countries  includes,  but  is  not  limited  to:  USA,  Germany,  Belarus,  Belgium,  Bulgaria,  Canada, 
Denmark,  Spain,  France,  Great  Britain,  Greece,  Hungary,  Iceland,  Italy,  Luxembourg,  Norway,  Netherlands, 
Poland,  Portugal,  Romania,  Russia,  the  Czech  and  Slovak  Republic,  Turkey,  Ukraine,  and  Georgia. 
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Recalling  the  commitments  they  have  made  in  the  Conference  on 
Security  and  Co-operation  in  Europe  to  promoting  greater  openness  and 
transparency  in  their  military  activities  and  to  enhancing  security  by  means  of 
confidence-  and  security-building  measures. 

Welcoming  the  historic  events  in  Europe  which  have  transformed  the 
security  situation  from  Vancouver  to  Vladivostok, 

Wishing  to  contribute  to  the  further  development  and  strengthening  of 
peace,  stability  and  co-operative  security  in  that  area  by  the  creation  of  an 
Open  Skies  regime  for  aerial  observation. 

Recognizing  the  potential  contribution  which  an  aerial  observation 
regime  of  this  type  could  make  to  security  and  stability  in  other  regions  as 
well. 

Noting  the  possibility  of  employing  such  a  regime  to  improve  openness 
and  transparency,  to  facilitate  the  monitoring  of  compliance  with  existing  or 
future  arms  control  agreements  and  to  strengthen  the  capacity  for  conflict 
preventions  and  crisis  management  in  the  framework  of  the  Conference  on 
Security  and  Co-operation  in  Europe  and  in  other  relevant  international 
institutions. 

Envisaging  the  possible  extension  of  the  Open  Skies  regime  into 
additional  fields,  such  as  the  protection  of  the  environment, 
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Seeking  to  establish  agreed  procedures  to  provide  for  aerial 
observation  of  all  the  territories  of  States  Parties,  with  the  intent  of  observing 
a  single  State  Party  or  groups  of  States  Parties,  on  the  basis  of  equity  and 
effectiveness  while  maintaining  flight  safety. 

Have  agreed  as  follows: 

1.  This  treaty  establishes  the  regime,  to  be  known  as  the  Open  Skies 
regime,  for  the  conduct  of  observation  flights  by  States  Parties  over  the 
territories  of  other  States  Parties,  and  sets  forth  the  rights  and  obligations  of 
the  States  Parties  relating  thereto. 

2.  Each  of  the  Annexes  and  their  related  Appendices  constitutes  an 
integral  part  of  this  Treaty  [Open  Skies  Treaty,  1992:  1-2]. 

1.2.  Background 

In  the  last  three  decades  most  academic  and  commercial  arenas  of  remote  sensing  have 
adopted  digital  data  processing  as  the  standard  paradigm  of  data  analysis.  Astronomers, 
geologists,  and  other  natural  scientists  are  trading  their  cameras  and  large  format  non-linear 
response^  photographic  plates  for  cameras  with  integrated  electronic  systems  composed  of 


^  Linear  and  non-linear  sensor  response  refer  to  the  relationship  between  the  input  and  the  output  of  the 
sensor.  Linear  sensors  have  outputs  which  have  a  constant  change  in  output  for  a  given  change  sensor  input 
over  the  design  range;  the  input  can  be  predicted  from  a  simple  linear  regression  model  as  a  function  of  the 
output.  Non-linear  sensors  cannot  be  described  as  a  simple  linear  functions  of  the  input  and  output.  The 
output  of  photographic  film  is  nonlinear.  Film  consists  of  an  acetate  strip  coated  with  silver  nitrate  grains 
boxmd  in  an  emulsion.  The  film  is  the  sensor  in  a  photographic  system,  where  the  input  is  the  irradiance  to 
which  it  was  exposed  for  a  fixed  amount  of  integration  time  and  the  output  is  the  emulsion  thickness.  The 
optical  density  (OD)  of  the  photographic  negative  increases  as  the  emulsion  thickness  increase.  This  OD  is 
used  as  the  measure  of  the  emulsion  thickness  and  it  is  proportional  to  the  logarithm  of  the  exposure  [Slater, 
1983:253-263]. 
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linear  opto-electronic  detectors,  such  as  charge  coupled  device  (CCD)  arrays  or 
photomultiplier  tubes,  and  high  powered  digital  computers  to  coordinate  the  data  collection 
and  analysis  for  both  infrared  (IR)  and  visible  imagery  data  [Menzel,  et  al,  1970:  94-96,  102], 
Data  stored  in  digital  format  can  be  transferred  across  network  or  telephone  communication 
lines,  reproduced  easily,  damage  resistant,  and  analyzed  by  multiple  parties  simultaneously. 
Conversely,  film  can  be  fragile  and  bulky,  only  reproduced  with  great  expense  or  loss  of 
information,  and  only  examined  by  one  group  of  analysts  at  a  time.  Even  with  the  advantages 
provided  by  digital  technology,  many  reasons  exist  for  analysts  to  continue  to  use 
photographic  film  for  aerial  imagery. 

Due  to  the  fall  of  the  Berlin  Wall  and  the  division  of  the  USSR  into  multiple  nation 
states,  many  countries  now  exist  whose  national  budgets  and  available  technologies  reinforce 
the  use  of  older  tools  to  perform  remote  sensing  tasks.  Aerial  cameras  using  high  resolution 
photographic  film  are  widely  available  technologies  with  which  aerial  reconnaissance  imagery 
data  has  been  collected  reliably  since  the  1940's.  One  notable  reason  for  this  film-based  view 
is  the  access  many  countries  have  to  the  photographic  reconnaissance  equipment  from  the 
Cold  War. 

Newer  European  sovereign  powers  need  the  same  strategic  information  about  their 
well  established  neighbors  to  enforce  their  nations’  security  as  their  neighbors  need  from  them. 
It  is  in  the  best  interest  of  the  more  technologically  advanced  countries  to  support  this 
gathering  of  information  through  remote  sensing  [AF  News  Service,  1995],  This  information 
increases  the  probability  of  long  lasting  stability  for  all  countries  by  increasing  the  confidence 
among  neighboring  states  that  no  hostile  movements  are  being  organized  {OST,  1992:  1-2], 
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As  stated  in  the  treaty  preface,  this  openness  was  the  central  intent  behind  the  Open  Skies 
effort,  and  lead  to  the  international  Open  Skies  Treaty.  Technological  equivalency  concerns 
led  the  signatory  countries,  or  State  Parties,  to  initially  accept  the  use  of  intelligence-type 
cameras  and  high  resolution  photographic  film,  not  digital  technology. 

The  OST  allows  signatory  countries  to  collect  aerial  imagery  of  strategic  areas  during 
flights  over  other  State  Parties’  borders.  For  the  ABM  Treaty  and  SALT  I,  verification 
imagery  was  limited  only  by  national  technical  means,  which  was  interpreted  by  the  US  and 
the  USSR  to  allow  the  use  of  any  means  of  verification  commanded  by  either  country 
[Luttwak,  1996].  The  US  and  the  USSR  possessed  evenly  matched  technologies  to  extract 
information  from  aerial  or  orbital  imagery.  The  level  of  technologies  commanded  by  members 
of  the  Open  Skies  regime  varies  to  such  an  extent  that  an  equitable  distribution  of  information 
could  not  be  guaranteed  without  administering  standard  restrictions  on  the  image  quality  each 
country  could  achieve.  Consequently,  to  implement  the  OST,  limitations  on  image  quality 
were  required.  The  difficulty  in  defining  an  actual  image  quality  metric  will  be  discussed  in 
Chapter  2,  but  as  the  matter  stands,  the  OST  contains  several  Decisions  defining  limits  on 
image  quality  as  factors  which  obstruct  the  ability  of  an  airborne  observer  to  see  ground 
objects  with  dimensions  smaller  than  30  cm  in  perfect  detail.  These  factors  are  synthesized 
into  a  function  of  the  flight  altitude,  the  length  of  the  smallest  identifiable  object  and  the  target 
modulation  [OST  Decision  14].  This  three  factor  combination  is  used  to  calculate  the 
minimum  altitude  at  which  an  OST  participating  aircraft  can  fly  over  an  area  of  interest  and  its 
result  is  computed  as  Although  the  identified  terms  in  have  not  yet  been  defined, 
for  reference  it  is  shown  in  Equation  (2.5),  it  reflects  the  value  placed  on  the  information 
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extractable  from  an  image  by  limiting  the  amount  of  ground  detail  any  one  participant  can 
achieve.  The  computation  of  H„i„  is  not  an  optical  performance  measure,  but  it  is  a  restriction 
on  data  collection  which  attempts  to  enforce  the  limits  placed  on  the  size  of  useful  details  in  an 
image. 

The  Open  Skies  Sensor  Working  Sub-Group  (OS  SWSG)  was  formed  as  the  technical 
arm  of  the  Open  Skies  (OS)  regime,  and  has  as  one  of  its  responsibilities  to  determine 
acceptable  methods  of  measuring  image  quality.  The  US  representatives  to  the  OS  SWSG 
report  directly  to  the  Office  of  Arms  Control  Implementation  and  Compliance,  Office  of  the 
Under  Secretary  of  Defense  for  Acquisition  and  Technology.  Along  with  representatives  from 
other  OST  countries,  they  chose  to  use  a  historically  sound  photogrammetric  approach  to  an 
image  quality  metric  composed  of  two  measured  factors. 

The  first  factor,  the  ground  resolved  distance  (GRD),  is  the  smallest  characteristic 
length  distinguishable  within  a  given  photographic  image  [Holst,  1995:216].  At  a  known 
altitude,  the  relationship  between  the  length  of  the  smallest  image  object  identified  on  a 
photographic  negative,  Axi,  and  length  of  the  actual  object  on  the  ground  are  related 
geometrically  as  shown  in  Equation  (1.1). 


where. 


GRD 


H,  is  the  altitude  or  distance  away  from  the  object  for  a  given  image, 
/  is  the  effective  focal  length  of  the  camera  system. 


(1.1) 
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GRD  is  a  common  measure  of  optical  system  performance  [Holst,  1995:  218],  The 
GRD  recorded  through  the  optical  system  increases  as  the  image  quality  degrades  (i.e.  with 
altitude  or  film  granularity).  The  GRD  for  a  given  OST  participating  airborne  optical  system 
is  estimated  by  photographing  a  ground  based  target  consisting  of  light  and  dark  bars  of 
decreasing  width  and  determining  the  width  of  the  smallest  discernible  bar  group.  The  value 
of  this  width  is  referred  to  as  L2  in  the  OST  documentation  [OST  Decision  14,  1994].  L2  only 
estimates  the  GRD  because  it  is  a  measurement  using  imperfect  detectors.  This  estimate  is 
only  an  upper  bound  on  the  GRD  because  the  bar  groups  occur  in  discrete  bar  widths.  If  the 
true  system  GRD  lies  between  two  of  these  sets  of  decreasing  width  bar  groups,  the  selection 
of  L2  will  always  be  greater  than  or  equal  to  the  actual  GRD.  This  measurement  technique 
will  be  discussed  in  depth  in  Chapter  2.  GRD,  however,  as  a  performance  measure  of  the 
optical  system,  is  not  a  measure  of  true  system  resolution,  but  a  measurement  of  the 
interaction  of  factors  such  as  the  system  optical  resolution,  atmospheric  effects,  object 
contrast,  and  sensor  sensitivity. 

The  second  component  to  the  OS  SWSG  image  quality  metric  is  the  target 
modulation.  Modulation  is  a  measure  of  the  range  in  exposure  values,  or  color  depth, 
contained  within  a  single  image.  All  details  in  an  image  are  described  on  the  sensor  as 
changes  in  exposure  [Holst,  1995:  216].  Exposure,  E,  is  a  measure  of  incident  energy  per  unit 
area'*  across  a  given  plane.  Modulation  can  be  expressed  empirically  in  terms  of  the  minimum 
and  maximum  exposure  values  within  an  image,  and  can  be  expressed  as  Equation  (1.2), 


'*  SI  units  of  E  are 


m 
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(E  -F  .  ) 

Modulation:  M  =  , 

(E  +E  ■) 

\  max  mm  / 


(1.2) 


where  Emax  and  Emm  are  the  maximum  and  minimum  values  of  exposure  within  an  image. 

When  using  photographic  film,  the  image  exposure  actually  recorded  on  film  is 
degraded  not  only  by  the  amount  of  irradiance^  lost  by  the  optical  train  in  diffraction, 
reflection,  refraction,  and  absorption,  but  also  by  the  ability  of  the  film  to  respond  to  small 
variations  in  the  exposure  distribution^.  The  film  response  and  the  irradiance  transmission 
through  the  optical  system  are  components  of  the  system  sensitivity.  Resolution  and 
sensitivity  are  separate  quantities.  Spatial  resolution'’  is  a  property  of  a  system  of  optics  and 
is  not  affected  by  noise,  film  response,  or  irradiance  transmission  through  the  system,  it  is  only 
a  function  of  the  geometry  of  the  optical  path  from  the  target  object  to  the  sensing  element 
(i.e.  the  film). 

The  measured  GRD  and  modulation  for  a  series  of  photographs  taken  over  a 
calibration  target  determines  the  altitude  an  Open  Skies  aircraft  can  fly  over  the  designated 
regions  of  interest  within  another  country.  The  OST  regime  has  chosen  to  use  ground  based 
targets  composed  of  painted  rectangular  panels  and  rectangular  bars.  An  Open  Skies 


^  Irradiance  is  the  rate  of  change  of  exposure  with  time  ( 2 ,  in  SI  units). 


m 


^  For  example,  if  the  smallest  visible  object  an  optical  system  can  resolve  has  length  10|a.m,  with  a  change  in 
exposure  of  1 ^ ,  but  the  film  can  only  record  changes  in  the  exposure  distribution  of  10  ,  the  GRD 


/m^ 


/m^ 


will  not  be  equal  to  the  optical  system  resolution,  because  the  object’s  contrast  is  not  great  enough. 

’  Three  other  aspects  of  resolution  can  be  described  [Holst,  1995:  218]; 

temporal  resolution  is  the  ability  to  separate  events  in  time— not  an  issue  with  single  frame  photography; 


gray  scale  resolution  is  a  measure  of  the  dynamic  range  (the  discernible  range  of  exposure  and  the 
smallest  detectable  change  in  exposure)  of  a  digitizer  for  data  storage  or  visualization  device  and  is 
strongly  affected  by  the  analog  to  digital  converter  design  and  noise  floor  of  the  electronics; 


spectral  resolution  is  the  range  of  the  electromagnetic  spectrum  discernible  by  a  system  (UV,  visible,  IR, 
etc.)  and  the  smallest  change  in  wavelength  it  can  detect. 
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calibration  target  can  have  many  configurations,  but  all  Open  Skies  calibration  targets  share 
the  series  of  light  and  dark-colored  bars  of  decreasing  bar  widths  with  separate  sets  in 
perpendicular  orientations. 

Figure  1  portrays  the  calibration  target  constructed  at  Wright-Patterson  Air  Force 
Base  (WPAFB)  for  the  US  Open  Skies  airborne  optical  system. 

1.2.1.  Calibration  Targets 

As  previously  discussed,  the  by-laws  of  the  OST  restrict  the  smallest  detail  in  an  image 
any  country  can  achieve  with  their  airborne  imagery  equipment  to  30  cm  for  an  acceptable 
target  modulation.  If  a  system  can  discern  an  object  with  a  GRD  less  than  30  cm  within  the 
accepted  modulation  range,  the  aircraft  must  fly  at  a  higher  altitude  to  reduce  the  ground 
resolution.  Each  optical  system  must  be  certified  in  each  of  its  configurations  to  be  used 
during  an  OST  flight.  The  image  quality  requirement  is  enforced  with  aerial  images  taken  of 
ground  based  calibration  targets.  These  targets  are  located  at  many  sites  around  the  world. 
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Figure  1,  Aerial  Photo  of  WPAFB  Calibration  Target  with  Annotated  Features 
(95-MOC-007 Mission,  Altitude  1697m,  KODAK  3404  Film) 
(Dimensions  represent  image  size  on  film  frame) 


An  aerial  image  of  the  target  constructed  at  Wright-Patterson  AFB,  OH,  from  OS 
Calibration  Mission  95-MOC-007  is  shown  in  Figure  1 .  The  important  features  of  the  target 
are  the  orthogonal  pairs  of  bar  triads  and  the  brightness  panels  and  are  indicated  in  Figure  1. 
Note  the  small  size  of  the  actual  image  digitized  from  the  photographic  film.  This  entire 
image  could  be  drawn  on  an  average  sized  fingernail.  The  small  size  of  the  image  explains  a 
majority  of  the  challenges  involved  with  the  development  of  an  automated  calibration  tool. 
Enough  samples  across  the  digitized  image  did  not  exist  to  use  common  image  identification 
tools  such  as  Fourier  analysis  or  matched  filter  techniques. 

The  paint  used  to  create  the  light  panel  and  the  light  bars  has  the  same  reflectance,  and 
a  lower  reflectance  paint  is  used  to  construct  the  dark  panel  and  dark  bars  with  equal 
reflectance.  The  reflected  light  energy  from  across  the  panels  is  used  to  measure  modulation. 
The  bars  are  constructed  with  decreasing  widths  and  separations  to  measure  the  GRD  of  aerial 
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photographic  systems.  These  types  of  ground  based  calibration  targets  degrade  over  time  due 
to  fading  of  the  paint  and  accumulation  of  environmental  pollutants.  This  degradation 
decreases  target  contrast  and  is  not  regularly  monitored  in  accordance  with  the  OST. 

The  largest  bar  group  set  has  a  width  of  120  cm  and  is  not  used  as  a  bar  group  for 
certification  purposes,  but  is  used  as  a  reference  width  in  this  thesis.  The  next  smaller  bar 
group  is  called  Bar  Group  #1  and  is  a  factor  of  (yfl  )^ ,  or  1 .26,  smaller  than  the  largest  bar 
group.  The  remaining  bar  groups  are  numbered  Bar  Group  #2  through  Bar  Group  #23  and 
decrease  in  width  by  a  V2  ,  or  1 . 12  ,  progression.  Bar  Group  #9  is  the  30  cm  wide  group 
and  is  marked  with  triangular  fiducials.  Both  orientations  of  the  bar  triads  are  used  to 
determine  the  smallest  GRD  and  the  brightness  panels  are  used  to  find  the  modulation.  The 
brightness  panels  and  bar  groups  are  both  required  to  have  a  modulation  between  0.66  and 
0.82  [OST,  Decision  14  1994:5]. 

1.3.  Problem  Definition 

Two  overwhelming  factors  motivate  the  development  of  automated  calibration 
techniques  for  the  OST  effort.  The  first  reason  supports  the  large  quantity  of  work  involved 
in  calibrating  each  sensor,  in  each  position  for  an  OST  mission.  At  any  time,  an  OST  country 
may  call  for  a  certification  flight  for  the  aircraft  mounted  optical  systems  [Open  Skies  Treaty, 
Decision  14].  Multiple  photographic  interpreters  (PI)  from  each  signatory  country  are 
responsible  for  determining  the  smallest  identifiable  bar  widths  by  viewing  dozens  of 
magnified  photographic  negatives  of  the  images  of  the  calibration  targets.  The  sheer  quantity 
of  images  which  must  be  used  in  order  to  obtain  the  values  defined  in  the  treaty  must  not  be 
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over  looked.  Ten  or  more  observers  from  participating  countries  could  be  employed  to 
critique  ten  or  more  bar  target  images  for  each  different  camera  system  to  be  flown.  Each 
country  is  permitted  several  dozen  flights  each  year  over  any  OST  nation.  The  sheer  quantity 
of  work  is  enough  to  justify  some  automated  simplification.  The  criteria  they  use  will  be 
discussed  in  Chapter  Two. 

The  second  motivating  factor  for  developing  an  automated  calibration  tool  involves 
the  reduction  of  uncertainty  found  in  many  common  calibration  techniques.  The  qualitative 
calibration  method  currently  employed  by  the  OST  signatories  is  influenced  by  human  factors 
such  as  experience,  fatigue,  and  ambient  light  level  [Slater,  1983:  241].  As  the  Open  Skies 
project  continues,  there  is  a  need  to  improve  consistency  in  the  data  processing  arena  ensuring 
no  one  country  gains  inequitable  advantages  over  others.  The  perception  of  an  unfair 
advantage  could  cause  conflict  among  OST  signatories  and  this  conflict  would  defeat  the 
purpose  of  the  treaty  and  its  philosophy  of  openness.  Quantitative  calibration  techniques 
would  reduce  uncertainty.  Hence,  the  creation  of  an  objective  measurement  could  increase 
the  probability  of  trust. 

Although  automated  computational  methods  of  measuring  imaging  system  capabilities 
are  being  created  in  Russia  and  other  countries,  personal  computer  power  and  digital  imagery 
are  not  widely  available  to  all  of  the  Open  Skies  participants,  and  as  such,  many  different 
methods  need  to  be  studied  and  debated  in  detail  to  increase  confidence  in  this  paradigm  shift. 
The  SWSG  has  been  directed  to  explore  automated  digital  methods  for  calibration  in  search  of 
one  which  can  be  readily  accepted  and  understood  by  all  Open  Skies  members. 
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This  thesis  supports  the  certification  automation  efforts  of  the  SWSG,  Ease  of  use  is 
required  for  any  viable  technique.  It  should  compute  valid  detection  parameters  given  only 
the  information  available  in  the  current  Open  Skies  data  collection  system.  All  algorithms 
must  run  on  widely  available  computer  equipment  and  use  commercially  available  software. 

Using  only  the  information  available  in  the  OSMPF  calibration  method  in  use  at  the 
time  of  this  study  reduces  the  number  of  analytical  procedures  available.  This  restrictive 
analysis  preference  constrains  the  problem  solution  to  techniques  which  only  depend  on  the 
knowledge  of  the  camera  geometry,  the  sensor  response,  the  flight  altitude,  the  physical 
dimensions  of  the  target,  and  the  characteristics  of  the  digitizing  scanner  employed.  No 
physical  transmission  model  of  the  optical  system  was  ereated  and  atmospheric  effects  were 
not  considered  directly.  Consequently,  no  knowledge  of  the  expected  signal  levels  is  available 
at  the  time  of  processing  is  available. 

The  central  issue  in  this  thesis  is  to  automatically  and  consistently  determine  which  bar 
triad  represents  the  smallest  set  resolved  for  all  photographic  images  of  the  calibration  target. 
More  than  one  method  to  satisfy  the  requirement  for  a  consistent  and  objective  digital  image 
processing  algorithm  was  pursued  in  this  thesis.  The  method  developed  determines  the  GRD, 
the  target  modulation,  and  the  minimum  altitude  at  which  an  OS  aircraft  can  fly,  These 
three  factors  were  calculated  and  compared  to  OSMPF  results  to  obtain  measures  of  merit  for 
the  algorithms  written  for  this  study. 
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In  the  method  used  in  this  thesis,  the  user  must  first  digitize  the  small  region  of  the 
photographic  film  negative^  containing  the  target  region.  Digitizing  the  photograph  creates  a 
very  large  bitmapped  file  which  contains  some  irrelevant  data.  Cropping  is  an  additional 
manual  step  which  is  useful  since  it  reduces  the  computation  time  of  the  analysis.  The  user 
then  applies  the  software  from  this  thesis  to  the  target  image  and  extracts  each  bar  group  from 
the  largest  to  the  smallest  widths  as  subimages.  The  software  processes  these  subimages  and 
collects  the  bar  group  statistics.  The  order  statistics^  of  the  digitized  image  gray  levels,  the 
gradients  of  the  data,  and  the  knowledge  of  the  target  layout  and  altitude  are  used  to 
determine  the  GRD.  The  modulation  is  determined  with  an  estimate  of  the  gray  level  values 
from  the  light  and  dark  brightness  panels,  which  are  converted  to  units  of  illuminance 
(Lumens)  using  calibrated  film  strips  of  known  exposure.  The  calibrated  film  strips  of  known 
exposure  are  referred  to  as  21  Step  Exposure  Wedges  which  are  film  strips  produced  at  the 
beginning  (header)  and  end  (tail)  of  each  roll  of  film  used  for  OST  flights.  The  21  Step 
Exposure  Wedges  are  exposed  with  known  exposure  using  Nation  Bureau  of  Standards 
calibrated  light  sources. 


*The  image  of  the  calibration  target  is  a  very  small  region  of  a  very  large  film  frame;  for  example,  in  the 
Kodak  3404  film  type,  the  frame  is  4.5”  square,  but  the  region  containing  the  calibration  target  is  less  than 
0.5”  square.  Digitizing  the  target  region  can  require  over  1  MB  of  8-bit  data,  so  even  with  many  megabytes  of 
RAM  on  a  computer,  it  is  best  to  extract  the  target  region  from  each  frame  first. 

^Order  statistics  are  used  in  this  thesis  because  the  distribution  of  the  bar  group  data  is  not  normal,  nor  are  the 
residual  values  from  the  estimated  model  of  the  bar  groups  used  to  calculate  bar  widths  [Kreyszig,  1983:  981], 
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1.4.  Summary  of  Key  Results 


Several  methods  for  evaluating  Open  Skies  imagery  were  developed  in  this  study  and 
tested  on  actual  data.  Due  to  an  undersampling  of  the  smaller  bar  groups  and  a  significant 
amount  of  noise,  Fourier  filtering  techniques  proved  to  be  unsatisfactory.  An  adaptive  digital 
technique,  however,  works  well  and  could  form  the  basis  for  automated  calibration  of  Open 
Skies  Treaty  systems.  This  method  matches  the  GRD  prediction  from  OSMPF  human 
observers  for  some  cases.  For  all  data  analyzed  in  this  thesis,  the  GRD  determined  by  ADiM 
was  smaller  than  or  equal  to  the  GRD  result  from  the  OSMPF.  This  interactive  technique  is 
adaptable  to  individual  bar  group  image  gray  level  amplitudes  and  is  based  on  the  order 
statistics  of  the  data  and  gradient  search  methods.  The  Adaptive  Digital  Method  (ADiM)  is 
implemented  within  a  graphical  user  interface  built  from  The  MathWorks  MatLab® 
programming  language  [MatLab®  Users  Guide].  Preliminary  results  suggest  a  reasonable 
match  with  human  observers.  Figure  2  and  Figure  3  demonstrate  the  correlation  between  the 
OSMPF  and  ADiM  computed  values  for 

Figure  2  shows  a  scatter  plot  of values  determined  by  the  OSMPF  (horizontal 
axis)  and  the  Hm„  values  determined  by  ADiM  (vertical  axis).  Figure  2  supports  the 
supposition  that  the  results  from  both  techniques  are  highly  correlated.  Creating  a  linear  fit 
between  the  two  sets  of  data,  containing  25  data  points,  in  the  form  of  a  simple  linear  model 
yields  a  correlation  coefficient  of  0.886.  This  means  almost  89  percent  of  the  variation  in  the 
OSMPF  Hmin  value  is  explained  by  the  variation  in  the  ADiM  value.  This  high  degree  of 
correlation  shows  that  the  compromises  made  in  the  adaptation  of  OSMPF  methods  to  digital 
data  can  be  used  to  produce  results  similar  to  the  manual  method. 
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for  OSMPF  vs  for  ADiM 


OSMPF  Hmin  (m) 

Figure  2,  Scatter  plot  o/Hmm  values  obtained  by  OSMPF  and  ADiM 


Explained  variation  does  add  to  the  credibility  of  the  results,  but  a  scatter  plot  of  both 
OSMPF  and  ADiM  computed  results  of//m,„  in  Figure  3  shows  that  none  of  the  ADiM  results 
exactly  matched  the  OSMPF  values  for  H„i„.  The  data  points  are  arranged  by  the  type  of  film 
used  for  data  collection.  Some  systematic  differences  in  the  two  methods  exist.  This 
difference  in  results  can  be  traced  most  often  to  different  values  of  the  calculated  film 
exposure  modulation  values  between  the  two  measurements. 
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Hmin  from  OSMPF  and  ADiM 


a 
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Data  Set  Number 

Figure  3,  Point  by  point  comparison  between  calculated  H min  values  for  OSMPF  and  ADiM 


A  reason  for  this  discrepancy  could  be  that  for  the  OSMPF  results  used  a  10  by  10 
sample  set  of  data  obtained  by  scanning  the  film  with  a  Perkin-Elmer  Digital  Scanner  (PDS)  to 
calculate  the  brightness  panel  modulation,  while  the  ADiM  method  used  data  digitized  with 
the  KODAK  5057  image  scanner  for  averages  of  hundreds  of  samples  all  over  the  brightness 
panels, 

A  simple  linear  model  was  built  with  the  least  squares  method  to  determine  the 
correlation  between  the  ADiM  and  OSMPF  Hmm  results.  In  order  to  assess  the  value  of  a 
simple  linear  model,  a  diagnostic  plot  of  the  residuals,  the  difference  between  the  estimated 
values  of  Hmm  from  the  model  and  the  actual  Hmm  values  calculated  by  the  OSMPF,  versus  the 
calculated  OSMPF  Hmm  values  shown  in  Figure  4  [Neter,  Wasserman,  Kunter  1990: 1 17],  The 
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abscissa  of  Figure  4  represents  the  difference  between  the  linear  model  and  the  OSMPF 
results.  The  three  different  symbol  demarcate  the  three  different  film  types  analyzed.  Note 
how  the  types  of  film  group  the  residuals  into  three  distinct  regions.  This  structure  could 
imply  that  higher  order  terms  are  needed  in  the  linear  model  to  explain  additional  variance  in 
the  OSMPF  results  with  the  ADiM  results. 


OSMPF  Hmin  Values  (m) 

Figure  4,  Linear  Model  Residuals  (Error)  plotted  against  OSMPF  Hmin 


Figure  4  demonstrates  the  possibility  that  different  films  affect  the  residuals  between 
measured  and  modeled  values.  If  it  were  important  to  convert  ADiM  values  to  OSMPF 
Hmin  values,  it  may  be  useful  to  include  a  variable  representing  film  type.  This  is  not  a 
requirement  and  will  not  be  discussed  further. 

This  structure  in  the  residual  plot  could  also  imply  that  improper  calibration  data  was 
used  to  compute  K2,  since  a  sample  of  each  type  of  film  must  be  used  in  order  to  compute  the 
modulation  and  translate  the  digitized  values  into  measures  of  luminance.  Further  study  is 
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required  to  fully  understand  the  contribution  to  the  errors  made  by  the  film  types.  However, 
the  correlation  of  88.6%  between  ADlM  and  OSMPF  does  support  that  the  software  methods 
developed  in  this  thesis  could  be  applied  to  the  OST  calibration  problem. 

1.5.  Organization  of  Thesis 

The  remainder  of  this  thesis  is  organized  as  follows.  Chapter  Two  reviews  the 
elements  of  optical  science  required  to  discuss  the  issues  of  image  quality  and  what  the 
measurements  of  the  OSMPF  actually  represent.  The  details  of  the  existing  methods  used  by 
the  OS  SWSG  to  determine  the  minimum  altitude,  //mm,  are  covered.  Chapter  Three  describes 
the  details  of  the  solution  method  and  the  Adaptive  Digital  Method  (ADrM)  algorithm 
developed  for  the  OSMPF  and  its  validity.  Tabulated  results  and  discussion  concerning  the 
results  are  contained  in  Chapter  Four,  followed  by  conclusions  and  suggestions  for  additional 
work  in  Chapter  Five.  MatLab®  code  listings  are  contained  in  appendices. 
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II.  Current  Methodology 

This  thesis  develops  a  computer-based  digital  calibration  method  to  determine  which 
bar  triad  has  the  smallest  discernible  width  in  a  digitized  photographic  imagery  of  an  OS 
calibration  target.  This  proposed  digital  method  is  derived  from  the  current  system  calibration 
method  employed  by  the  signatory  countries  to  the  OST,  which  was  developed  from  the  need 
to  assure  all  OS  participants  received  the  same  detail  of  information.  They  have  agreed  to  a 
level  of  image  quality  based  on  two  quantities;  the  ground  resolved  distance  (GRD)  and  the 
target  modulation,  defined  in  Equation  (1.2)  [OS  Decision  14,  1994].  The  OS  calibration 
technique  and  their  image  quality  criteria  are  discussed  in  the  following  pages  and  the  details 
of  the  digital  method,  ADiM,  will  be  the  focus  of  the  next  chapter. 

The  OST  and  the  appended  Decisions  do  not  refer  to  the  selected  criteria  as  measures 
of  image  quality.  In  spite  of  this  omission,  the  two  factors  identified  by  the  OST  used  as  a 
standard  for  image  information  content  regulation  actually  do  constitute  a  true  image  quality 
metric  and  restrict  the  information  recorded  in  an  image  by  defining  the  lowest  altitude  at 
which  an  OST  participant  can  record  image  data  over  another  OST  country.  Much  more  will 
be  discussed  about  this  in  Section  2.2.1.  Image  quality  is  a  judgment  on  the  inherent  value  of 
the  information  content  in  an  image,  and  H„i„  restricts  the  level  of  detail  recordable  by  an 
airborne  remote  sensing  system.  The  process  of  this  information  extraction  is  called  image 
analysis,  whose  principles  encompass  the  methods  for  “detecting,  identifying,  and  measuring 
objects  of  interest  from  the  aerial  perspective,”  according  to  the  Manual  of  Remote  Sensing 
[Estes,  et  al.,  1983:987]. 
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The  first  factor  in  the  OST  image  quality  equation  is  the  smallest  width  discernible  in 
an  aerial  image  taken  of  a  ground  target,  at  a  known  altitude,  and  it  is  represented  by  a 
measure  of  the  minimum  GRD,  called  Z,^.  The  second  factor  is  the  film  incident  image 
exposure  modulation,  called  K2.  Although  the  OST  has  no  formal  definition  of  an  image 
quality  metric,  in  using  these  factors  to  limit  the  flight  altitude  of  the  data  collection  aircraft, 
the  OS  regime  has  implicitly  created  a  metric.  However,  the  GRD  is  not  an  independent 
component  of  image  quality,  but  a  function  of  the  target  contrast  as  well  as  the  width  of  the 
target  on  the  ground  and  the  geometric  resolution  of  the  camera  optics. 

GRD  performance  is  measured  for  each  aerial  imaging  system  in  use  by  the  Open 
Skies  participants,  in  each  of  the  mounting  locations  and  configurations  to  be  certified  for  a 
data  collection  flight.  GRD  is  called  resolution  in  the  OST;  however,  in  this  thesis,  the 
ground  width  of  the  smallest  set  of  bars  deemed  distinguishable  will  be  called  the  GRD.  The 
term  resolution  can  refer  to  four  different  physical  performance  properties  of  an  image 
acquisition  system  as  discussed  in  Section  2.2.  The  measurement  of  GRD  is  affected  by  many 
environmental  factors  and  the  strongest  parameters  are  illustrated  in  Figure  5.  The  GRD 
increases  as  the  light  reflecting  off  of  the  passive  target  is  transmitted  through  the  atmosphere. 
This  increase  is  caused  by  the  loss  of  light  due  to  an  unfavorable  incident  light  angle  with  the 
sun  as  the  day  passes  and  the  seasons  change,  image  smearing  due  to  imperfect  camera  shutter 
speed  timing,  imperfections  in  the  bar  target  surface  paint,  atmospheric  temperature  changes 
which  distort  the  reflected  image  as  it  propagates  through  to  the  sensor,  and  finally  sensor 
imperfections.  Smearing  of  the  images  in  the  direction  of  motion  of  the  aircraft  also  occurs 
due  to  finite  exposure  time  of  the  film. 
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Figure  5,  Sources  of  GRD  Degradation^^ 

The  terms  distinguishable  and  resolved  are  interpreted  in  the  OST  to  mean  the 
smallest  bar  group  on  which  both  dark  bars  are  differentiable  from  their  adjacent  light  bars, 
along  their  entire  length.  Figure  6  portrays  the  details  of  each  bar  to  be  distinguished.  By 
definition  in  the  OST,  a  bar  group  is  defined  as  a  pair  of  light  bars  separated  by  a  dark  bar  of 
the  same  size.  Hence,  in  Figure  6  the  three  light  bars  depicted  represent  two  triple  bar  groups 
where  the  left  most  light  bar,  the  center  light  bar,  and  the  dark  area  between  them  represent 
the  left  most  triple  bar  group.  The  dark  area  between  them  is  considered  the  dark  bar.  The 
center  bar  and  the  right  most  bar  group  define  the  outside  bars  of  the  right  most  triple  bar 
group.  The  bar  length  of  both  light  and  dark  bars  is  five  times  longer  than  the  bar  width.  The 
use  of  this  target  will  be  discussed  in  more  detail  later  in  this  chapter.  This  standard  for 
distinguishability  will  be  relaxed  for  the  digital  method  discussed  in  Chapter  Three,  because 
true  edges  in  the  image  are  difficult  to  define  in  the  presence  of  noise  and  the  digitized  data 
has  additional  noise  that  could  make  a  group  of  bars  appear  indistinguishable  when  the 
information  content  of  each  bar  is  actually  resolvable. 


^'’The  GRD  of  the  bar  triad  reflected  irradiance  can  never  increase  as  it  propagates  along  the  optical  path  to 
the  receiving  sensor. 
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Figure  6,  Triple  Bar  Group  Diagram 


GRD  for  a  photographic  camera  system  is  a  function  of  the  resolution  of  the  system, 
the  response  of  the  film  to  light,  the  irradiance  incident  on  the  ground  from  the  sun,  the 
particulate  content  of  the  air,  and  the  modulation  of  the  calibration  target,  to  name  a  few 
[Holst,  1995:  241],  Hence,  using  the  GRD  as  part  of  the  image  quality  criteria  considers  the 
interaction  of  the  of  all  of  these  quantities  for  the  given  data  run,  not  any  of  the  individual 
factors  directly.  In  order  to  further  discuss  the  strengths  and  weaknesses  of  the  OST 
approach,  the  concepts  introduced  above,  resolution,  sensitivity,  and  modulation  must  be 
established. 

2.1.  The  Basic  Physics  of  OST  Image  Quality 

For  all  studies  which  require  an  accurate  and  reproducible  measurement  of  image 
quality,  the  actual  definition  of  an  image  quality  criterion  is  difficult  [Holst,  1995:  239]. 
Defining  the  overall  goal  of  the  data  collection  should  drive  the  details  of  this  metric.  A  large 
variety  of  measures  exist  to  describe  the  information  content  of  a  recorded  image,  and  for  the 
most  part  these  criteria  are  not  interchangeable  among  different  sensors  and  may  not 
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correspond  to  engineering  or  scientific  measures  of  spatial  resolution,  dynamic  range,  or 
spectral  resolution  [Holst,  1995:  239],  For  instance,  the  metric  for  observing  a  tan, 
camouflaged  truck  in  the  tan,  brush-covered  desert  under  a  Saudi  Arabian  summer  sun  at  a 
distance  of  two  kilometers,  is  not  the  same  system  performance  needed  for  seeing  a  black 
truck  on  the  Siberian  tundra  driving  across  the  snow  white  tundra  from  an  airplane  at  noon. 

The  important  distinction  between  the  two  scenarios  is  that  the  former  is  mostly  an 
issue  of  detector  sensitivity  and  the  latter  is  mostly  one  of  system  resolution.  It  is  also 
important  to  recognize,  as  will  be  discussed,  the  effects  of  sensitivity  and  resolution  are 
complimentary.  They  are  easily  confused  in  determining  image  quality  because  the  eye  has 
difficulty  distinguishing  between  a  loss  of  sensitivity  or  a  loss  of  resolution.  Both  effects  can 
cause  a  loss  in  image  detail  to  a  detector  or  a  human  observer. 

Before  it  is  an  image  recorded  on  some  type  of  receiver,  such  as  film  or  electronic 
detectors,  a  light  source  has  created  a  signal  which  has  been  reflected  off  or  transmitted 
through  several  different  types  of  media.  The  ground,  the  atmosphere,  lenses,  and  mirrors  can 
absorb  and  distort  portions  of  this  image  forming  light  reducing  the  strength  and  warping  the 
geometry  of  the  original  signal.  The  film  or  electronic  detectors  have  ranges  of  signal  strength 
sensitivity  over  which  they  can  react  and  record  incident  light  to  form  an  image  recognizable 
to  a  human  observer.  If  the  incident  signal  is  too  weak  or  the  optical  path  has  absorbed  too 
much  of  the  signal,  the  detector  can  not  react,  and  therefore  it  can  not  form  the  details  of  the 
image.  If  the  signal  is  too  strong,  the  detector  becomes  saturated  and  can  no  longer 
distinguish  between  incremental  changes  in  the  signal  strength  which  would  form  the  details  of 
an  image.  These  are  the  issues  of  sensitivity.  Furthermore,  not  only  the  optics  of  an  optical 
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system  attenuate  available  signal  strength.  Film  and  electronic  detectors  also  have  signal 
degrading  natures  and  act  as  attenuating  and  distortion  filters  as  the  transmitted  image  passes 
onto  them.  Film  has  image  recording  capabilities  which  depend  on  the  size  of  the  smallest 
grain  of  emulsion  and  on  the  smallest  amount  of  light  it  takes  to  change  the  color  intensity. 
Other  factors  affect  the  performance  of  film  and  Slater  discusses  them  in  detail,  but  they  will 
not  be  discussed  in  this  thesis.  The  film  used  for  the  OST  can  have  a  sample  size  near  1  nm, 
which  represents  a  ground  size  of  approximately  1  mm  at  a  flight  altitude  of  2000m.  The 
KODAK  5057  used  in  this  thesis  can  have  sample  sizes  near  bpm. 

The  effective  optical  system  aperture  also  acts  as  a  detail  reducing  filter.  These 
degrading  filter  effects  further  complicate  the  attempt  to  define  an  image  quality  parameter  by 
the  fact  that  all  images  can  be  considered  to  be  composed  of  an  infinite  number  of  frequencies, 
but  the  optical  system  can  only  capture  a  small  fraction  of  these  frequencies.  Sharp  edges  are 
composed  of  nearly  infinite  numbers  of  the  highest  frequencies,  so  for  bar  target  based 
methods,  the  first  details  to  be  lost  are  the  sharp  edges  along  the  perimeter  of  each  bar  in  a  bar 
group.  In  order  gain  a  greater  understanding  of  the  effect  an  optical  system  can  have  on  a 
reflected  image,  the  system  modulation  transfer  function  must  be  discussed. 

2. 1. 1.  The  system  modulation  transfer  function. 

To  some  extent,  even  physics  works  against  the  definition  of  a  universal  image  quality 
metric.  This  especially  affects  passive* '  techniques  using  bar  targets  or  other  edge  defining 
measurement  methods.  Not  only  do  the  bars  of  higher  spatial  frequency  have  less  surface  area 
to  reflect  incident  light  than  their  neighboring  bars  of  larger  widths  (and  lengths  for  the 


*  *  Passive  imaging  techniques  have  no  source  of  lighting  beyond  the  natmal  light  provided  by  the  sun. 
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WPAFB  target),  but  the  finite  diameter  optics  further  attenuate  the  narrow  bars  due  to 
diffraction,  resulting  in  the  attenuation  of  higher  spatial  frequencies  described  by  overall 
system  spatial  modulation  transfer  function  (MTF)  [Castleman,  1995:371-377],  TheMTF 
measures  the  attenuation  of  spatial  frequencies  due  to  the  physical  geometry  of  the  optics,  but 
does  not  show  any  information  about  the  sensitivity  of  the  system  or  its  detector  [Holst,  1995: 
237], 

The  spatial  variations  in  a  recorded  image  depend  on  the  geometry  of  the  object(s)  in 
the  image  and  the  geometry  of  the  data  collection  optical  system,  The  shape  and  size  of  the 
optical  system  as  a  whole  determines  the  smallest  distance  resolvable  in  the  exposure 
distribution  incident  on  the  detector.  The  study  of  Fourier  Optics  [Goodman,  1988 
(Reissued)]  treats  the  subject  of  electro-magnetic  radiation  propagation  through  optical 
systems  in  depth.  In  his  classic  text,  Goodman  develops  the  basic  concept  of  the  MTF  and  the 
relationship  between  the  object  spectrum,  0(f),  the  image  spectrum,  1(f),  and  MTF  given  by: 

I{f)  =  MTF^,Uf)0{f)  (2,1) 

The  system  MTF  is  a  useful  performance  measure.  It  represents  the  attenuation  of  spatial 
frequencies,  which  compose  the  details  visible  in  an  object,  due  to  the  system  being  analyzed. 

The  MTF  is  measured  using  the  image  a  point  source  of  light  Avith  known  signal  strength. 

Because  each  element  in  along  the  optical  path  of  an  image  contributes  to  the  over  all  system 
MTF  and  can  be  modeled  as  series  of  linear  systems  such  that  the  effect  on  the  image  by  the 
atmosphere,  the  lens  system,  and  the  windows  and  mirrors  encountered  by  the  propagating 
image  can  be  represented  in  the  following  fashion: 
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where,  ,  is  the  transfer  function  describing  the  attenuation  of  detail  due  to  the 

atmosphere,  ,  the  transfer  function  for  the  aircraft  window  where  the  image 

enters  the  optical  system,  ,  the  transfer  fimction  for  the  optical  path  along  which 

the  image  transmits,  and,  MTF^i^(f) ,  the  transfer  function  for  the  emulsion  on  the  film. 

The  MTF  of  a  data  collection  system  is  a  filter  fianction  which  determines  how  spectral 
components  of  the  object  are  attenuated  or  boosted  as  they  pass  through  the  particular  system 
of  interest.  If  the  MTF  is  thoroughly  known,  then  all  of  the  original  object  details  could  be 
completely  reconstructed  from  the  recorded  image  by  dividing  the  MTF  from  both  sides  of 
Equation  (2.1)  and  retrieving  the  details  of  the  original  object.  Unfortunately,  the  OSMPF 
does  not  have  the  end-to-end  MTF  curves  for  their  optical  systems,  and  the  atmospheric  MTF 
would  need  to  be  measured  for  each  different  image  in  order  to  remove  its  effects. 

Optical  systems  are  only  one  category  of  systems  for  which  a  filter  function  describes 
the  effect  the  system  has  on  signals  in  the  frequency  domain.  The  human  ear  is  also  a 
sophisticated  data  collection  system.  Sound  waves  are  attenuated  by  the  human  ear  such  that 
most  humans  carmot  hear  waves  with  frequencies  higher  than  20  kHz  or  lower  than  100  Hz. 
For  comparison  of  an  ear  to  an  optical  system,  such  as  a  camera,  different  spatial  frequencies 
are  manifested  in  the  sizes  of  the  objects  to  be  imaged.  Higher  spatial  frequencies  appear  as 
small  objects  or  fine  detail  in  an  image,  lower  spatial  frequencies  correspond  to  large  objects. 

The  on-axis  MTF  for  the  US  Open  Skies  KS-87  lens  was  measured  by  the  Wright 
Laboratory’s  Metrology  Lab  and  is  close  to  the  diffraction  limit.  This  lens  was  used  in  the 
collection  of  all  of  the  data  analyzed  in  this  thesis.  However,  no  theoretical  analysis  of  the 
total  system  MTF  was  performed  due  to  a  lack  of  information  about  the  MTF  of  the  film 
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digitizer.  The  lack  of  measured  local  atmospheric  conditions  also  prevented  an  accurate 
theoretical  analysis. 

For  most  imaging  tasks,  the  primary  goal  of  collection  is  to  gain  or  record  useable 
information.  No  direct  correlation  exists  between  the  information  derived  from  an  image  and 
the  shape  of  the  acquisition  system’s  modulation  transfer  function  (MTF)  [Holst,  1995:  239; 
Slater,  1983:  240-241].  Nonetheless,  human  appreciation  of  an  image  increases  when  the 
signal  irradiance  level  is  high,  the  fine  details  in  the  image  exist,  and  the  image  graininess  is 
low.  Humans  tend  to  judge  an  image  of  this  nature  to  have  better  image  quality  than  one  of 
lower  signal  strength  or  higher  noise  [Holst,  1995:  239].  Systems  with  high  MTF  values  at 
high  spatial  frequencies  are  capable  of  providing  highly  detailed  images. 

The  MTF  of  the  optical  system  does  not  hold  the  answer  to  the  question  of  finding  a 
single  number  to  describe  an  image  quality  metric.  It  is  only  a  tool  to  describe  the  effect  of 
the  optical  system  on  the  final  image.  Since,  the  exposure  distribution  of  the  resulting 
recorded  image  is  reduced  by  the  same  proportional  amount  at  each  spatial  frequency,  the 
final  value  of  an  image  is  still  a  function  of  the  original  details  it  contained  and  the 
interpretation  of  the  user. 

2. 1.2.  Modulation  and  Contrast  Ratio  Are  Different  Photogrammetric 
Values. 

The  terms  contrast  and  modulation  are  well  defined  photogrammetric  quantities  and 
both  quantitatively  compare  measurements  of  the  maximum  and  minimum  exposure  recorded 
on  the  film  or  other  sensor.  For  the  OST  image  analysis,  the  averaged  digitized  gray  level 
values  of  the  brightness  panels,  shown  in  Figure  5,  are  used  to  compute  the  maximum  and 
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minimum  logarithm  of  exposure  used  to  compute  the  contrast  and  modulation.  Contrast  is 
shown  in  Equation  (2.3)  and  modulation  is  defined  in  Equation  (2.4)  [Slater,  1983:241]. 


Contrast  Ratio:  C  = 


(2.3) 


where,  Emax  is  the  maximum  exposure  value  and  is  the  minimum  exposure  value. 
Modulation  is  a  measurement  of  the  difference  between  the  maximum  and  minimum  exposure 
level  normalized  by  the  sum  of  the  two  and  contrast  is  simply  the  ratio  between  the  two 
exposure  extremes.  The  algebraic  relationship  between  modulation  and  contrast  is  seen  in 


(2.4). 


Modulation.  M  = - 

C  +  1 

The  value  of  modulation  is  bounded  between  zero  and  one.  However,  the  contrast 
ratio  is  an  unbounded  quantity  which  only  describes  the  ratio  of  exposures,  the  original 
exposure  level  is  masked. 

Figure  7  demonstrates  the  results  of  contrast  and  modulation  for  a  signal  with  a 


(2.4) 


maximum  and  minimum  difference  of  500-^  at  five  different  signal  levels.  Notice  that  a 

m 

measure  of  exposure  near  zero  causes  the  contrast  ratio  value  to  increase  without  bound,  but 
the  modulation  simply  approaches  unity.  The  graph  shows  that  the  Modulation  values  are 
bounded  between  zero  and  unity,  yet  the  Contrast  Ratio  is  unbounded  as  the  difference 
between  the  minimum  and  maximum  exposure  increases. 


30 


Table  1,  Demonstration  of  Contrast  and  Modulation  Values  for  a  Fixed  Exposure  Di  fference 


E„ax(^) 

E„i„  {tL) 

Contrast  Ratio 

Modulation 

500 

1 

500 

0.9960 

1000 

500 

2,0000 

0.3333 

10,000 

9500 

1.0526 

0.0256 

100,000 

99,500 

1.0050 

0.0025 

1,000,000 

999,500 

1.001 

0.0003 

Demoastration  of  Boimdedness  of  Modulation 


Figure  7,  Comparison  of  Modulation  and  Contrast  for  a  Fixed  Exposure  Difference  of 500 pJ 


The  modulation  calculated  within  an  image  and  the  actual  modulation  value  in  the 
MTF  are  not  the  same  physical  quantities.  An  important  distinction  between  M  in  Equation 
(2.4)  and  the  modulation  value  in  the  MTF  is  that  M  is  defined  from  exposure  differences 
within  a  single  image,  but  the  MTF  values  are  defined  as  comparisons  of  exposure  as  a 
function  of  spatial  frequency  at  the  input  and  output  of  an  optical  system.  Therefore,  in  this 
thesis  the  actual  MTF  cannot  be  calculated  from  the  available  data  because  no  input  irradiance 
measurements,  across  any  band  of  spatial  frequencies,  were  available  for  comparison. 
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2. 1. 3.  Under  sampling  reduces  the  certainty  of  identifying  individual  bar 
groups 

The  undersampling  in  the  digital  image  scanner  which  hinders  the  bar  width 
measurement  can  be  observed  in  Table  2.  The  sample  size  for  the  data  set  from  the  95-MOC- 
007  mission  was  30.98  microradians/pixel  wide,  which  projects  to  8.88  cm/pixel  on  the 
ground  at  an  altitude  of  2750m,  using  KODAK  3412  film.  For  bar  groups  8,  9,  10,  and  11, 
Table  2  shows  that  the  width  of  each  bar  ranges  from  2.7  to  3.8  samples  wide.  Distinguishing 
between  them  with  only  the  estimated  pixel  widths  could  be  difficult  since  they  each  are 
sampled  with  at  most  three  full  pixels  for  this  sampling  resolution.  As  discussed  later, 
averaging  improves  the  estimate  of  the  bar  widths  and  permits  ADiM  to  effectively  distinguish 
between  bar  groups  8  through  1 1 .  Bar  groups  8  and  9  were  not  as  differentiable  as  the  other 
groups  and  this  can  be  explained  by  noticing  that  the  bar  widths  of  bars  8  and  9  are  only 
different  by  0.4  pixels,  both  only  3  full  pixels  wide.  Bar  group  8  has  a  ground  width  of  3.8 
and  bar  group  9  is  3.4  pixels  wide,  even  using  averaging,  the  variance  of  the  data  was  too 
large  to  consistently  differ  between  them. 

The  uncertainty  due  to  undersampling  does  not  affect  the  accuracy  or  usefulness  of 
ADiM-it  just  demonstrates  that  the  OS  Treaty  criteria  applied  to  the  data  acquisition  system 
in  use  at  the  writing  of  this  thesis  cannot  be  easily  fully  automated.  Man  in  the  loop  analysis 
was  required.  Notice  how  small  the  bar  group  widths  are  on  the  film  image.  However,  if  the 
criteria  and  data  acquisition  were  redesigned  in  order  to  eliminate  the  man  in  the  loop,  total 
computer  automation  could  be  possible.  Potential  designs  are  discussed  in  Chapters  4  and  5. 
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Table  2,  Comparative  Number  of  Samples  Across  Bar  Groups  from  One  Target  to 
Demonstrate  Source  of  Uncertainty  as  the  Number  of  Samples  Across  Each  Bar  Group 

Decrease 


Bar  Group 

Width  on  Film  (cm) 

Angular  Subtense  @  2750m 

Samples 

1 

0.004877 

75.6 

275  (32.3) 

8.5 

2 

0.00436 

67.4 

245 

7.6 

3 

0.003873 

60.0 

218 

6.75 

4 

0.003442 

53.5 

194.6 

6.0 

5 

0.003098 

47.6 

173 

5.4 

6 

0.002725 

42.4 

154 

4.75 

7 

0.002433 

37.8 

137 

4.24 

8 

0.00218 

33.7 

122.6 

3.8 

9 

0.001951 

30 

109 

3.4 

10 

0.001721 

26.7 

97.1 

3.0 

11 

0.001549 

23.8 

86.6 

2.7 

2.2.  Current  Procedures 

The  Open  Skies  Media  Processing  Facility,  OSMPF,  operated  by  the  US  Air  Force 
National  Air  Intelligence  Center  (NAIC)  at  WPAFB,  has  among  its  many  duties  the  task  of 
data  management  for  data  taken  from  the  US  Open  Skies  aircraft  and  each  of  its  sensors.  To 
certify  a  single  sensor  in  a  particular  configuration,  they  compile  the  PI  recommended 
minimum  resolvable  bar  widths,  referred  to  as  L2,  from  several  aerial  passes  over  the  ground- 
based  calibration  targets  and  calculate  the  minimum  height  at  which  an  OS  aircraft  can  fly  with 
that  sensor.  They  also  work  wdth  other  country’s  OS  offices.  This  effort  results  in  many 
redundant  tasks,  hundreds  of  rolls  of  film,  and  megabytes  of  data  to  manipulate  and  store. 

2.2.1.  Calibration  Target  Data  Collection  and  Analysis 

Regardless  of  the  calibration  target  configuration,  the  OS  aircraft  image  acquisition 
system  to  be  tested  flies  over  the  target,  its  optical  systems  record  the  target,  the  resulting 
images  are  processed,  and  the  actual  acetate  film  negatives  are  given  to  the  PI  for  evaluation 
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on  a  light  table  under  a  microscope.  The  PI  then  determines  which  sets  of  bars  are  the 
smallest  resolvable  in  each  flight  direction.  This  minimum  bar  width  information  along  with 
the  exposure  modulation  of  the  brightness  panels  is  used  to  determine  the  minimum  altitude 
the  particular  aircraft  can  fly  for  its  image  gathering  mission  discussed  in  later  sections. 
Determination  of  the  smallest  bar  groups  and  flight  altitude  is  done  by  vote  among  the  PI. 

The  smallest  bar  width  belongs  to  the  bar  group  for  which  80  percent  of  the  PI  involved  voted 
to  have  the  smallest  resolved  bar  width.  The  value  of  L2  in  the  equation  is  assigned  the 
width  meeting  this  80  percent  rule  criterion. 

The  treaty  does  not  require  a  single  standard  configuration  for  a  calibration  target,  but 
size  and  reflectance  restrictions  do  exist.  These  specifications  are  discussed  in  section  1.2.1. 
According  to  the  treaty,  these  targets  will  all  consist  of  pairs  of  light-dark-light  tri-bars,  also 
referred  to  as  bar  triads,  or  light-dark  bi-bars  as  discussed  earlier  and  shown  in  Figure  6. 

The  pair  of  bar  triads  or  the  bi-bars  will  be  referred  to  as  a  bar  group.  The  methods  in  this 
thesis  only  were  tested  on  bar  triad  type  targets  and  would  require  software  updates  to  be 
used  on  other  targets.  With  the  proper  digitization,  there  is  no  intrinsic  reason  the  recording 
system  must  be  an  optical  camera  and  film  system.  For  example,  the  method  developed  in  this 
thesis  is  based  on  the  geometry  of  the  bar  target  and  adapts  to  the  statistical  nature  of  each 
target  as  will  be  discussed  in  Chapter  3.  The  sets  of  three  light  bars  shown  in  Figure  6 
actually  comprise  a  pair  of  bar  triads  as  the  treaty  specifications  are  stated.  The  two 
orientations  of  the  bar  triads  record  spatial  frequency  in  directions  along  or  in  (parallel  to)  the 
aircraft  line  of  flight  and  across  (perpendicular  to)  the  aircraft  line  of  flight. 


The  data  analyzed  in  this  thesis  is  based  on  a  target  whose  bar  triad  aspect  ratio 
remains  a  constant  of  5: 1  (length  to  width).  Figure  1  depicts  an  aerial  photo  of  the  actual 
target  used,  which  is  one  of  the  USA’s  OS  calibration  targets,  and  is  located  on  Wright- 
Patterson  Air  Force  Base,  in  Ohio.  In  addition  to  the  bar  groups,  the  targets  also  include 
brightness  panels  of  much  larger  dimension  than  even  the  largest  bar  triad,  and  have  the  same 
reflectance  properties  as  the  bars  and  will  be  discussed  in  detail  later. 

The  participating  PI  view  an  image  frame  through  a  microscope  on  a  light  table.  Both 
dark  bars  from  each  triad  must  be  visible  to  deem  a  set  of  bars  resolved.  The  decision  process 
to  select  the  accepted  smallest  bar  set  for  a  given  pass  leads  to  a  collation  of  the  chosen  bar 
groups  from  each  PI  and  then  the  choice  is  made  based  on  the  smallest  set  of  bars  on  which  at 
least  80  percent  of  the  PI  agree.  Each  PI  votes  on  their  preferred  bar  group  set  and  the 
corresponding  GRD.  The  smallest  GRD  with  80  percent  of  the  total  votes  determines  the 
final  choice  of  resolved  length.  This  80  Percent  Rule  generally  yields  reasonable  results  for 
Hmin,  but  is  not  an  analytical  measure  of  optical  system  resolution  and  hence  it  is  subjective. 

The  minimum  altitude  an  aircraft  can  fly  with  a  given  optical  system,  Hmm,  is  found 
with  two  pieces  of  information  from  the  recorded  ground-based  target:  1)  the  smallest  set  of 
bar  triads  that  can  be  visually  separated,  and  2)  the  exposure  contrast  between  the  light  and 
dark  brightness  panels  found  on  the  target.  The  brightness  panels  of  the  target  have  the  same 
reflectance  as  the  bright  and  dark  portions  of  the  bar  triads,  and  their  average  modulation  is 
currently  measured  using  a  single  machine  which  combines  a  Zeiss  microdensitometer  with  a 
PDS  (the  machine  is  referred  to  as  the  Micro-D)  to  digitize  the  brightness  panels'^.  Two 

A  KODAK  5057  high  resolution  photographic  scanner  is  currently  being  tested  as  an  alternative  to  the 
PDS,  but  as  of  Jun  ’96,  the  PDS  was  still  the  primary  source  of  brightness  panel  illuminance  measurement. 
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calibration  steps  are  required  to  go  from  the  output  of  the  PDS,  or  any  scanner,  to  find  the 
measurements  in  the  terms  needed  to  process  the  data. 

As  the  first  step,  the  PDS  is  used  to  digitize  21 -Step  wedges  on  processed  film 
exposed  by  National  Bureau  of  Standards  (MBS)  calibrated  light  sources  to  calibrate  the  PDS 
output  values  of  analog-to-digital  units  (ADU)  to  optical  density  (OD).  The  PDS  is  operated 
in  a  mode  such  that  the  output  ADU  are  linearly  proportional  to  OD.  A  curve  can  then  be 
interpolated  resulting  in  a  direct  linear  conversion  of  ADU  to  OD  for  the  PDS.  To  confirm 
the  film  was  exposed  in  valid  exposure  ranges  and  to  do  the  altitude  calculation,  the  next 
calibration  step  will  convert  the  data  into  units  of  the  light  exposure  incident  on  the  film — so 
one  more  step  is  required.  OSMPF  has  the  capability  to  generate  calibrated  strips  of  film  on 
which  the  exposure,  realty  the  logarithm  of  the  exposure  (log£ )  is  known  to  high  accuracy. 
These  strips,  called  21 -step  exposure  wedges,  are  created  on  the  beginning  (head)  and  end 
(tail)  of  each  roll  of  film  processed. 

The  wedges  are  then  scanned  with  the  PDS,  the  head  and  tail  results  are  averaged,  and 
the  resulting  ADU  versus  \ogE  curve  is  transformed  to  a  sigmoid  shaped  OD  to  logii  .  The 
resulting  curve  is  called  the  Hurter-Driffield  (H-D)  Curve,  or  the  D  log  £  curve.  A  sample 
curve  generated  for  the  NBS  and  OSMPF  wedges  is  shown  below  in  Figure  8. 
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Hurter-Driffield  Curve  for  Kodak  3412  Film  Emulsion 


Log(Exposure) 

Figure  8,  Hurter-Driffield  Curve  For  KODAK  3034  Film,  with  fitted  line  to  the  linear  region 

The  digitized  portions  of  the  brightness  panels  are  then  transformed  into  units  of 
logE  and  those  points  are  plotted  on  the  H-D  curve  to  make  sure  they  fall  in  the  linear  region 
of  the  curve.  This  hnear  region  is  demarcated  by  the  least  squares  fit  straight  line  in  Figure  8 
using  “x”  symbols  to  mark  the  predicted  optical  density.  Currently,  the  points  used  to 
determine  this  linear  fit  are  determined  by  the  analyst. 

As  was  mentioned  in  Chapter  One,  a  different  digitization  method  was  used. 
Specifically,  the  data  was  scanned  with  a  KODAK  5057  photographic  film  scanner.  At  the 
time  of  the  writing  of  this  document,  the  OST  participants  had  not  accepted  the  validity  of  the 
use  of  the  KODAK  5057  scanner  in  the  calibration  mission  procedures.  The  only  difference  in 
the  two  methods  is  the  type  of  scanner  used.  The  critical  difference  between  the  two 
instruments  is  that  the  PDS  is  an  NBS  accredited  instrument  accepted  widely  in  the 
photometric  community,  and  its  linearity  is  calibrated  with  NBS  instruments. 
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The  image  data  is  then  used  to  calculate  Hmm  with  a  Microsoft  Excel®  and 
VisualBASIC®  macro/program  called  K2.xls.  K2.xls,  is  used  with  the  PI  recommended 
minimum  bar  width,  labeled  L2 ,  in  the  Hmi„  formula  [OST,  Decision  14]: 


n 


i=l 


'30' 

'0.40‘ 

U2J 

0.45 


(2.5) 


E 


(2.6) 


C-1 


C  +  1 


(2.7) 


where, 

n  ,  number  of  passes  made  over  the  target, 

H.  ,  the  height  above  ground  level  at  which  the  data  was  collected,  in  meters 

Zj  ,  the  ground  separation  of  two  adjacent  bars  in  the  minimum  resolvable  bar  triad,  in 

centimeters 

E  ,  the  exposure  represented  by  the  average  density  in  emulsion  over  a  brightness  panel, 
AlogE  ,  the  difference  of  the  logE  value  between  the  light  and  dark  brightness  panels. 
C  ,  is  the  measurement  of  panel  contrast,  units  of  exposure. 

The  expression  for  K2  in  Equation  (2.7)  is  the  OST  label  for  the  physical  value  of 
modulation,  M,  expressed  previously  in  Equation  (2.4).  Hence,  the  modulation,  M,  will  be 
referred  to  as  K2. 


2.2.2.  Bar  Targets  Intrinsically  Add  Uncertainty  to  GRD  Calculation 

The  OSMPF  approach  to  determining  L2  is  based  on  historically  accepted  methods  for 
judging  the  quality  of  an  aerial  photographic  camera  system  in  the  lab,  under  collimated,  fully 
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illuminated  conditions  [Slater,  pg  241],  Slater  also  notes  that  the  original  technique  specified 
the  use  of  a  transparent  slide  of  the  bar  target  in  transmission  from  the  focal  plane,  not  from 
reflectance  as  in  the  case  of  the  ground  based  targets.  He  warns  that  even  in  these  controlled 
conditions. 

The  subjective  and  statistical  nature  of  resolving  power  should  not  be 
overlooked,  nor  should  the  fact  that  the  yfl  increments  in  the  size  of  the 
standard  elements  in  a  target  correspond  to  a  change  of  12  percent  in  the 
spatial  frequency.  The  reader  may  often  be  uncertain  over  a  range  of  two  or 
three  target  elements. 

Hence  using  this  criterion  on  data  collected  under  uncontrolled  yield  GRD 
determinations  which  are  larger  than  the  GRD  which  more  accurately  describes  what  level  of 
detail  an  imaging  system  can  record,  Human  factors  such  as  fatigue  or  political  agendas  could 
fiirther  complicate  this  issue. 

The  selection  of  a  bar  group  with  width  larger  than  the  true  system  GRD 
underestimates  the  system’s  true  performance  and  will  produce  an  altitude  of  at  least  12 
percent  lower  than  it  should  be.  This  is  due  to  the  discrete  %l2  bar  width  progression 
multiplier  applied  to  the  bar  widths  in  the  construction  of  OS  bar  targets  and  also  due  to  the 
fact  that  Hmirt  is  linear  in  Once  an  automated  digital  technique  is  introduced  as  an 
acceptable  data  process,  errors  introduced  in  the  digitization  process  will  further  complicate 
GRD  determination  by  human  observers.  Degradation  in  the  digitized  image  is  called 
quantization  error''*,  and  the  image  will  be  further  degraded  due  to  noise  introduced  via  the 
electronics  required  for  digitization. 


(Note:  %l2  =  1.12246) 

Quantization  error  or  noise  is  a  result  of  taking  a  continuous  signal  and  breaking  it  into  discrete  blocks  of 
values.  It  is  a  form  of  aliasing  which  is  using  data  taken  with  a  sample  size  larger  than  the  smallest  detail  in 
the  original. 
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A  key  issue  for  both  human  and  digital  methods  of  GRD  estimation  is  the 
determination  of  where  the  light  and  dark  bars  are  separated.  The  OSMPF  defines  the  two 
adjacent  light  and  dark  bars  to  be  separated  if  there  exists  notable  light  and  dark  regions  along 
the  entire  length  two  bars.  Unfortunately,  using  a  digital  method  this  criteria  is  not  adequate 
due  to  the  blurring  and  noise  corruption  of  boundaries.  This  blurring  systemically  occurs 
from  the  finite  lens  diameter  and  noise  is  introduced  by  the  electronics  and  by  the  physics  of 
the  light  detection  process  [Roggemann  &  Welsh,  1996].  However,  in  averaging  within  the 
same  image,  much  of  the  random  noise  is  attenuated  leaving  a  picture  of  the  bar  group  cross 
section  which  is  easier  to  interpret.  This  averaging  can  be  along  the  length  of  the  bar  group  of 
a  single  image,  as  done  in  ADiM,  or  across  multiple  images  of  the  target  at  the  same  altitude 
and  within  the  same  region  of  airspace.  So,  averaging  can  also  increase  the  details  available 
within  an  actual  data  run  if  enough  frames  in  the  same  airspace  are  taken. 

The  next  nine  image  are  presented  to  demonstrate  some  of  the  imaging  process  affects 
on  the  bar  group  images  as  bar  widths  decrease.  This  presentation  also  alludes  to  the  subject 
matter  of  the  next  chapter,  which  discusses  the  development  of  ADiM  and  the  idiosyncrasies 
in  the  data  which  make  ADiM  challenging.  The  following  bar  group  data  are  fi-om  OS 
Calibration  Mission  095-MOC-007  (Pass  1,  Frame  16).  All  of  the  following  data  was 
digitized  from  this  mission  taken  over  the  WPAFB  target  at  approximately  1700  meters  above 
the  target.  The  following  figures  depict  three  different  bar  groups  in  three  different  views. 
The  presentation  of  these  three  sets  of  figures  further  demonstrates  the  ways  in  which  a 
human  observer  could  be  convinced  that  information  was  not  available  within  an  image,  when 
in  fact,  the  important  detail  was  in  the  photograph  buried  beneath  the  image  noise. 
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The  first  figure  in  each  set  will  be  the  gray  scale  image  representation.  The  second 
figure  displays  the  image  as  a  surface  plot  so  the  smoothness  of  the  edges  of  each  bar  group 
can  be  see,  and  the  third  figure  of  each  set  contains  two  representations  of  the  lateral  mean 
cross  section  of  each  bar  group  and  a  10*^  order  polynomial  fit  of  that  cross  sectional  plot. 
Figure  9  defines  the  average  of  the  lateral  mean  cross  section  of  a  bar  group  image.  This  is 
another  way  to  represent  digital  imagery  data  and  takes  advantage  of  the  statistical  error 
reduction  attributed  to  arithmetic  averaging. 


Figure  9,  Demonstration  of  Averaging  Direction  to  Calculate  the  Lateral  Mean  Cross- 

Section  of  a  Bar  Group  Image 

The  imaging  effects  on  the  bar  groups  are  caused  by  the  optical  system  MTF  and 
losses  due  to  noise  or  other  random  effects.  A  few  of  these  properties  can  be  clearly  observed 
in  the  following  figures.  Among  these  important  properties  are:  the  blurring  of  the  bar  edges, 
the  diminishing  peak  bar  group  amplitude,  and  the  assimilation  of  the  dark  bars  by  the  spread 
of  the  light  bars.  These  are  the  important  effects  of  the  system  MTF  on  the  images  in 
providing  digital  tools  to  the  OSMPF.  Figures  10,  1 1,  and  12  highlight  the  details  of  bar 
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group  #2  (barwidth=67.35  cm-ground^^)  from  this  data  set.  The  succeeding  figures  present 
bar  groups  9  and  10. 

Figure  10  is  a  gray  scale  image  of  bar  group  #2  and  shows  the  light  and  dark  bars  with 
clearly  defined  edges,  although  those  edges  are  not  sharp.  Each  pixel  represents  a  6  pm 
sample  from  the  digitized  film  frame.  The  spread  of  the  light  bars  already  encompass  much  of 
the  dark  bar  area  notable  in  this  image.  Bar  group  #2  is  the  second  largest  bar  group  and  the 
MTF  is  already  making  the  task  of  edge  detection  difficult.  Consequently,  difficulty  in  edge 
detection  translates  into  uncertainty  in  determination. 
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Figure  1 1  is  a  surface  plot  of  bar  group  #2  and  confirms  any  suspicions  that  this  bar 
group  image  is  indeed  that  of  a  pair  of  resolved  bar  triads.  The  distinction  between  the  light 
and  the  dark  bars  has  been  eroded,  although  it  does  still  exist.  The  strength  of  the  light  bar 
signal  is  demonstrated  for  bar  group  #2.  This  surface  plot  of  bar  group  #2  also  shows  the 
noise  apparent  in  the  data.  The  cause  of  the  jagged  peaks  is  not  explained  by  this  image,  but 


5  io  15  20  25  30  35 

Figure  10,  Grayscale  Image  Bar  Group  #2 


Barvvidth  as  measured  on  the  ground  target. 
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causes  could  stem  from  the  atmospheric  scintillation,  imperfections  in  the  paint  surface  across 
the  bars,  or  any  of  the  other  systems  involved  in  making  this  image. 
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Figure  11,  Surface  Plot  of  Bar  Group  #2 

The  next  representation  of  bar  group  #2,  in  Figure  12,  consists  of  two  plots,  a  solid 
line  connecting  data  points  composing  the  average  lateral  cross-section  of  bar  group  #2  and 
the  dashed  line  tracing  the  10*^  order  polynomial  fit  [MatLab®  Reference  Manual]  found  for 
these  data  points.  Again,  two  distinct  pair  of  bar  triads  are  apparent.  The  asterisks  present  on 
the  central  peak  represent  the  centroid  of  the  cross  section  and  the  two  points  adjacent  to  it. 
These  markers  were  used  in  ADiM  to  initiate  a  peak  and  valley  identification  search  which  will 
be  discussed  in  Chapter  3 . 
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P1F16.TIF.  Group  #2  .  Across  LOF 


Figure  12,  Mean  Lateral  Cross-Section  Plot  of  Bar  Group  #2  and  lO"'  Order  Polynomial  Fit 
(Solid  Line  :  Bar  Group  Data  ;  Dashed  Line:  lO'^  Order  Polynomial  Fit;  Markers:  the 

centroid  and  adjacent  points) 


In  the  digitized  image  and  surface  plot  of  bar  group  #9  (Barwidth=30.00  cm-ground) 
in  Figures  13  and  14,  one  cannot  have  complete  confidence  in  distinguishing  each  edge  along 
any  bar  edge  in  the  gray  scale  image.  The  noise  masks  the  information  in  the  data. 
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Figure  13,  Grayscale  Image  Bar  Group  #9 


Figure  14,  Surface  Plot  of  Bar  Group  #9 


However,  by  averaging  along  the  rows,  as  shown  in  its  mean  lateral  cross  section  in 
Figure  15;  the  light  and  dark  bars  are  clearly  distinguished  and  the  polynomial  fit  further 
clarifies  the  distinction  between  the  two  bar  triads.  In  the  preceding  Figure  1 1,  bar  group  #2, 


the  dark  bar  amplitude  values  appear  to  be  near  the  amplitude  values  of  the  background 
surrounding  the  bar  group  surface.  But  for  bar  group  #9,  the  dark  bar  amplitudes  are 
significantly  higher  than  the  peak  amplitudes  from  the  light  bars.  As  the  bar  group  widths 
decrease  in  the  image,  the  light  bars  continue  to  spread  until  the  dark  bar  amplitudes  become 
difficult  to  discern.  In  Figure  15,  the  mean  lateral  cross  section  of  bar  group  #9,  it  is  evident 
that  the  dark  bar  amplitudes  are  significantly  above  the  back  ground  amplitude  level.  Figure  9 
demonstrates  how  this  cross  section  is  generated.  This  is  an  effect  of  the  imaging  system 
MTF  on  the  reflected  bar  image  and  will  be  discussed  in  detail  in  the  next  section. 


P1F16  TIF,  Group  #9  ,  Across  LOF 


Figure  15,  Mean  Lateral  Cross-Section  Plot  of  Bar  Group  #9  and  lO'^  Order  Polynomial  Fit 
(Solid  Line:  Bar  Group  Data  ;  Dashed  Line:  Order  Polynomial  Fit;  Markers:  the 

centroid  and  adjacent  points) 

For  the  next  smaller  bar  group,  bar  group  #10  (barwidth=26  cm-ground),  no 
definitive  bar  edges  can  be  detected  in  the  gray  scale  image  shown  in  Figure  15,  and  again,  no 
edges  can  be  detected  in  the  surface  plot  in  Figure  17. 
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Figure  16,  Grayscale  Image  Bar  Group  #10 


Figure  17,  Surface  Plot  of  Bar  Group  #10 


As  seen  in  the  images  for  bar  group  #9,  with  the  use  of  the  noise  reduction  properties 
from  averaging,  clear  separations  between  light  and  dark  bars  are  easily  distinguishable  and 


are  demonstrated  in  Figure  18,  despite  the  fact  that  this  information  was  not  available  in  either 
of  the  previous  two  plots  for  bar  group  #10. 


P1F16.TIF,  Group  #10  ,  Across  LOF 


Figure  18,  Mean  Lateral  Cross-Section  Plot  of  Bar  Group  U 10  &  id’'  Order  Polynomial  Fit 
(Solid  Line:  Bar  Group  Data  ;  Dashed  Line:  Kf’'  Order  Polynomial  Fit  Markers:  the 

centroid  and  adjacent  points) 

2.2.3.  Optics  distort  GRD  measurement,  and  surface  roughness  distorts 
target  modulation  measurement 

Not  only  is  the  determination  of  the  L2  measurement  of  the  GRD  estimation  adversely 
affected  by  measurement  noise  and  target  inconsistencies,  the  target  modulation  value,  K2, 
measurement  is  distorted  as  well.  OSMPF  computes  the  illuminance  of  the  light  and  dark 
brightness  panels  by  averaging  a  10  pixel  by  10  pixel  region  manually  extracted  from  each 
panel.  One  hundred  samples  are  produced  and  the  average  pixel  value  from  each  panel  is 
converted  to  Emax  and  and  used  to  compute  K2. 
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Figure  19  shows  a  30  pixel  by  60  pixel  subsection  sample  of  the  light  and  dark 

brightness  panels  and  shows  the  inconsistency  of  the  reflected  surface  and  suggests  the  need 

for  a  larger  sample  set  to  calculate  the  actual  modulation  of  a  target. 
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Figure  19,  Surface  Plot  of  Brightness  Panels  for  OS  Mission  095-MOC-007  to  Emphasize 

Surface  Inconsistencies 

Its  light  panel  values  range  from  230  ADU  to  160  ADU  and  its  dark  panel  values 
range  from  100  ADU  to  150  ADU.  The  following  histogram  depicts  the  actual  spread  of  the 
data  values,  for  the  entire  brightness  panel  region. 
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Brighness  Panels  Histogram 


Figure  20,  Histogram  of  Brightness  Panel  to  Show  Value  Variance 
Figures  19  and  20  suggest  the  need  for  a  larger  sample  size  than  10  pixels  by  10  pixels, 
currently  in  use  by  the  OSMPF,  In  developing  the  ADiM  software,  fluctuations  in  the  K2 
value  as  much  as  25  percent  were  noted  by  simply  changing  the  size  and  location  of  the 
sample  set  across  the  brightness  panels. 

2.3.  Chapter  Summary 

In  order  to  construct  the  OST  and  gain  the  consensus  of  each  of  the  desired 
participants,  limits  were  placed  on  the  smallest  details  recordable  by  airborne  optical  systems. 
These  limits  defined  acceptable  sizes,  L2,  and  modulations,  K2,  of  details  within  an  image  and 
were  used  to  calculate  the  minimum  height  an  aircraft  could  fly  to  collect  imagery  under  the 
hospices  of  the  OST.  Effectively,  the  OST  participants  defined  quantities  whose  combination 
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was  used  to  determine  a  measure  of  image  quality.  Image  quality  defines  the  value  of  the 
information  contained  within  an  image. 

Image  quality  definitions  are  inherently  subjective  due  to  the  complicated  interactions 
between  human  perception,  image  transmission,  and  image  detection.  The  OST  image  quality 
standard  set  limits  on  the  smallest  detail  recordable  and  the  largest  range  in  exposure  levels 
observable  in  a  image  of  a  bar  target.  Unfortunately,  the  smallest  width  observable  within  an 
image  of  the  original  ground  target  is  a  function  of  both  sensitivity  and  resolution. 

Additionally,  the  modulation,  or  exposure  range,  within  an  image  is  also  a  function  of  both  the 
sensitivity  and  resolution  of  the  entire  optical  system. 

For  a  real  image,  the  sensitivity  factor  is  a  limitation  if  the  irradiance  level  of  the 
reflected  light  from  the  target  transmitted  through  the  optical  path  of  the  airborne  system  is 
not  strong  enough  or  the  film  is  not  sensitive  enough  to  react  to  and  record  the  exposure  from 
the  incident  image.  The  resolution  factor  usually  describes  the  attenuation  rate  of  an  optical 
system  on  incident  image  spatial  frequency  content.  If  the  signal  strength  of  the  smallest 
spatial  frequencies  transmitted  from  the  ground  object  and  passing  through  the  airborne 
optical  system  to  the  film  is  not  strong  enough  to  be  detected  by  the  film,  no  record  of  those 
details  will  be  observable  in  the  final  image.  In  other  words,  if  the  MTF  of  the  entire  optical 
system  has  attenuated  spatial  frequencies  available  to  the  system,  but  reduced  their  exposure 
below  the  sensitivity  of  the  film,  the  details  described  by  those  spatial  frequencies  will  not  be 
observable  in  the  final  image. 

Therefore,  the  two  parameters  used  to  compute  Hmm  are  not  independent.  In  fact, 
with  a  thorough  characterization  of  the  optical  system  and  atmospheric  measurements  for  each 
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image,  the  MTF  for  the  total  system  could  be  reproduced  and  removed  from  the  digitized 
image.  This  new  image  reproduces  the  ground  object  almost  perfectly.  This  could  cause 
future  problems  for  countries  whose  technology  does  not  allow  them  the  ability  to  make  these 
post-processing  corrections  and  would  deny  them  details  that  more  technologically  advanced 
countries  could  record  from  them. 

The  correlated  OST  criteria,  the  discontinuous  construction  of  the  bar  targets,  and  the 
coarseness  of  the  digitization  of  the  bar  target  images  lead  to  a  measure  of  image  quality 
which  is  difficult  to  automate.  The  inconsistent  surfaces  of  the  bar  target  and  the  unknown 
background  noise  level  make  the  development  of  a  fully  automated  bar  group  finder  very 
difficult,  because  the  image  on  the  film  frame  is  very  small  and  the  orientation  of  the  bar  target 
is  not  consistent  from  frame  to  frame.  The  inability  to  find  edges  around  the  smaller  bars  due 
to  noise,  the  physical  rounding  effects  of  the  imaging  process,  and  the  insufficient  numbers  of 
samples  across  bar  widths  around  30  cm-ground  and  smaller  makes  the  implementation  of  the 
OST  bar  resolution  standard  impossible  to  implement. 

Even  with  a  perfectly  automated  method,  with  the  current  technique  of  data 
acquisition  and  digitization,  limitations  on  the  confidence  in  the  measurement  of  the  GRD  and 
modulation  are  unavoidable.  The  fact  that  the  bar  targets  are  constructed  with  decreasing 
discrete  bar  widths  reduces  confidence  in  the  GRD  measurement,  because  the  actual  value  of 
the  GRD  will  almost  always  exist  between  two  bar  groups.  The  analyst  is  only  able  to  say  the 
system  GRD  is  no  larger  than  a  given  bar  group  width,  but  can  never  bound  the  GRD  result 
with  tolerances  smaller  than  the  difference  between  two  sets  of  bar  groups. 
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The  digital  method  developed  in  this  thesis  uses  a  set  of  criteria  which  is  based  on  the 
spirit  of  the  OST  standards  of  image  quality,  but  does  not  replicate  their  technique, 
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III.  The  Adaptive  Digital  Method  (ADiM) 

As  discussed  previously,  the  optical  system  performance  criteria  used  by  the  OS 
regime  depends  on  two  components:  the  target  illuminance  modulation,  K2,  and  the  target 
GRD,  L2.  These  parameters  determine  the  minimum  altitude,  Hmm  ,  at  which  an  Open  Skies 
aircraft  can  fly.  K2  is  measured  directly  from  the  brightness  panels  on  the  target  in  all 
calibration  methods.  Any  value  of  K2  greater  than  0.4  is  penalized  by  increasing  this  can 
be  observed  in  the  K2  definition  in  the  previous  chapter.  L2  is  measured  subjectively  according 
to  criteria  concerning  the  geometry  of  the  bar  triads,  and  the  OS  regime  agreed  value  for  Z2  is 
30  cm  [OST,  Decision  14].  Chapter  2  contains  the  discussion  of  the  component 
calculation  details. 

This  tool,  the  Adaptive  Digital  Method  (ADiM),  allows  the  user  to  extract  the  target 
image  features  required  to  estimate  signal  strength,  contrast,  and  location  against  the  target 
background  environment.  The  accuracy  of  these  estimates  will  determine  the  validity  of  the 
calculated  values  for  L2  and  K2.  Because  these  limits  are  placed  on  the  values  for  K2  and  L2 
to  accomplish  the  goals  of  the  OST  by  its  participants,  the  inferred  quantity,  Hmi„,  represents  a 
value  judgment  on  the  information  content  of  the  image.  Hence,  Hmi„  is  a  measure  of  image 
quality.  Holst  has  taken  great  pains  to  emphasize  the  subjective  nature  of  image  quality 
[Holst,  pg  239].  In  fact,  with  calibration  targets  such  as  the  WPAFB  target,  all  subjectivity  in 
the  OS  certification  process  cannot  be  removed  [Slater,  pg  246].  The  intent  of  ADiM  is  to 
offer  the  OS  another  tool  with  which  to  manage  large  quantities  of  data  and  increase  the 
certainty  held  for  the  H„i„  computation. 
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3.1.  Technique  Overview 


ADiM  applies  standard  digital  processing  techniques  to  digitized  files  created  fi'om 
photographic  negatives  to  determine  the  orientation  and  amplitude  modulation  of  individual 
bar  groups.  First  and  second  order  difference  functions  are  used  in  ADiM  for  gradient  search 
techniques  and  edge  detection.  The  flow  chart  in  Figure  21  describes  the  steps  in  the  ADrM 
process.  Bar  target  analysis  with  ADiM  is  iterative  and  requires  human  interaction.  No 
consideration  is  made  in  the  current  version  of  ADiM  for  rotational  transformations,  so  the 
user  must  ensure  that  the  rows  and  columns  of  the  digital  image  lay  parallel  or  perpendicular 
to  the  orientations  of  the  bar  groups  of  interest. 

ADiM  operation  is  simple  and  can  be  summarized  in  four  easy  steps:  setup,  select, 
analyze,  and  repeat  as  needed.  The  user  selects  a  bar  group,  in  either  orientation.  If  the 
software  can  determine  the  orientation  of  the  selection,  and  guess  the  bar  group  number,  with 
its  corresponding  width,  of  the  selection  within  one  or  two  bar  widths,  the  bar  group  is 
considered  resolved.  ADiM  is  then  inactive  until  the  next  bar  group  is  selected.  This  routine 
is  followed  until  the  software  can  no  longer  recognize  the  orientation  of  the  bar  group.  In 
development,  it  was  observed  that  if  the  orientation  was  not  detected  correctly,  then  the  bar 
width  was  also  reported  incorrectly. 
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3.1.1.  Step  by  Step  ADiM 

A  data  run  begins  in  ADiM  the  same  way  it  begins  in  the  OSMPF  technique.  Both 
ADiM  and  the  OSMPF  digitize  the  21 -Step  Exposure  Wedges,  also  known  as  sensitometric 
strips,  and  store  the  digitizer  output,  in  ADU  values,  to  a  computer  file.  The  output  files  of 
the  sensitometric  strip  data  are  regressed  against  their  known  Exposure  values,  using  a  simple 
linear  model.  From  the  results  of  this  computation,  the  conversion  factors  from  ADU  values 
to  units  of  Exposure  are  estimated,  and  used  without  error  analysis.  A  more  detailed 
description  of  this  calibration  and  the  use  of  photographic  film  is  discussed  later  in  this 
chapter.  The  conversion  equation  coefficients  are  recorded  and  this  completes  the  block 
diagram  step,  Use  Sensitometric  Strip  to  Calibrate  Gray  Level  Curve.  ADiM  then  requires 
the  digitization  of  the  bar  target  region  under  analysis.  The  region  containing  the  bar  target  is 
scanned  with  the  smallest  sample  size  available  on  the  digitizer. 

Before  launching  the  ADiM  software,  the  user  must  update  the  MatLab^  routines  with 
the  results  of  the  gray  level  to  logarithm  of  exposure  calibration,  the  flight  altitude,  the 
effective  focal  length  of  the  optical  system,  and  the  sampling  size  at  which  the  data  was 
digitized.  The  user  must  also  enter  the  name  used  to  identify  the  data  run.  This  completes  the 
Update  Program  with  Pass  Parameters  step  of  the  block  diagram. 

ADiM  assumes  the  target  images  are  available  as  digitized  image  files  in  tagged  image 
file  format  (TEFF).  In  the  current  version  of  the  software,  these  parameters  must  be  entered 
directly  into  the  MatLab^^  code.  In  future  versions,  this  information  should  be  placed  directly 
in  the  header  of  the  TIFF  image  as  part  of  the  data  storage  process.  This  procedure  completes 
the  Setup  Complete,  Load  Digital  Image  stage  of  the  ADiM  block  diagram. 
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ADiM  does  not  model  the  physics  or  optics  of  the  image  acquisition  system  nor  does  it 
take  into  account  the  effects  of  image  transmission  through  the  atmosphere.  This  technical 
data  is  not  readily  available  and  would  require  alterations  in  the  current  OST  data  acquisition 
technique.  For  that  reason,  the  solution  presented  here  uses  only  knowledge  of  the  target 
geometry  and  the  stochastic  properties  of  the  digitized  image  data,  which  is  readily  available 
from  the  OSMPF. 

The  user  now  starts  MatLab®  and  then  launches  the  ADiM  macros.  The  ensuing 
description  follows  the  Determine  L2  branch  of  the  block  diagram.  ADiM  enables  the  user  to 
determine  a  target  GRD,  L2,  in  each  of  the  two  perpendicular  bar  group  orientations.  The 
GRD  is  defined  as  the  width  of  the  smallest  bar  group  for  which  the  software  can  identify  the 
bar  group  orientation^^  and  calculate  a  valid  estimate  of  the  bar  group  width  and  geometry. 
After  the  ADiM  routines  have  been  launched,  two  dialog  boxes  appear  enabling  the  use  to  use 
a  previous  file  or  select  a  new  one.  If  a  new  file  is  being  examined,  a  standard  MS  Windows 
dialog  box  appears  to  enable  the  user  to  browse  the  computer  file  system  to  find  the  correct 
file.  The  user  selects  the  filename  of  the  bar  target  to  be  analyzed.  If  the  file  selected  is  new 
or  was  previously  analyzed,  the  user  sees  the  screen  shown  in  Figure  22. 


This  criterion  is  the  same  as  that  used  by  Professor  Beloglazov  in  his  1993  discussion  of  applying  least 
squares  indicators  to  the  OS  problem.  The  methods  to  determine  the  orientation  and  bar  width  in  this  thesis 
are  unique  to  the  best  knowledge  of  the  author  [Beloglazov,  1993]. 
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Figure  22,  Opening  Screen  to  ADiM 

In  the  opening  screen  there  are  five  important  features.  The  most  prominent  feature  is 
the  image  of  the  entire  bar  target  along  the  left  side  of  the  window.  The  remaining  features 
consist  of  the  upper  menu  bar  across  the  top  of  the  window  and  the  three  control  buttons  on 
the  lower  left  half  of  the  widow.  The  upper  menu  bar  allows  the  user  to  exit  and  print  the 
current  window.  The  use  of  the  control  button  icon  is  necessary  to  engage  the  processing 
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algorithms  of  ADiM.  The  ZOOM  button  allows  the  user  to  enlarge  regions  of  the  bar  target 
image  using  point  and  drag  mouse  actions  to  make  region  selection  easier.  After  the  area  of 
interest  is  in  view,  the  ZOOM  button  is  pressed  again  to  disable  the  zoom  feature  (the  button 
toggles  the  zoom  function  on  and  off).  The  SELECT  REGION  h\}\Xon  allows  the  user  to 
select  the  region  of  interest  (ROI)  containing  the  bar  group  and  to  launch  the  ADiM  routines 
for  analysis.  For  a  typical  data  iteration,  the  user  would  use  the  ZOOM  control  to  enlarge  the 
area  containing  the  ROI,  toggle  off  the  zoom  fimction,  and  press  the  SELECT  REGION 
control  to  select  the  ROI. 

ADiM  then  captures  the  selected  region  and  presents  it  to  the  user  for  verification 
shown  in  Figure  23 .  A  pop-up  window  appears  and  queries  the  user  as  to  the  status  of  the 
presented  selection.  If  the  image  is  the  intended  ROI,  the  user  selects  OK,  the  processing 
begins,  and  the  results  are  shown.  If  the  image  is  not  the  intended  ROI,  he/she  rejects  the 
query,  and  the  selection  process  starts  again.  This  is  an  iterative  process,  especially  for  the 
smaller  bar  groups. 

3. 1.1.1.  Finding  L2 

The  process  to  determine  L2  is  started  by  selecting  the  Select  Region  option  and 
dragging  a  rubber-band  box  around  the  region  containing  the  bar  group  or  brightness  panel  to 
be  analyzed.  Figure  23  displays  the  result  of  the  selection  of  bar  group  #1,  for  the  25  Aug  95 
mission. 
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Figure  23,  Selected  Bar  Group  #1  Using  ADiM 


This  ends  the  block  diagram  step.  Select  Region  Containing  Bar  Group. 

The  next  two  steps  along  the  block  diagram  denote  progress  toward  extracting  the  bar 
group  of  interest  from  the  ROI.  These  steps  are  labeled:  Median  Filter  Bar  Group  Subimage 
and  Isolate  Bar  Group  from  Background.  The  initial  ROI  is  surrounded  by  unnecessary 
background  data  which  adds  uncertainty  to  the  parameters  to  be  estimated.  It  is  useful  to 
isolate  the  bar  group  pixels  from  those  containing  background  data.  The  group  perimeter,  or 
a  defined  boundary  around  the  outside  of  the  bar  group  must  be  found  to  isolate  the  bar  group 
from  the  background.  Two  MatLab^  IPT functions  are  used  to  find  this  perimeter,  bwimg 
and  perimeter. 

To  find  this  group  perimeter,  a  binary  image  is  created.  A  binary  image  is  an  image 
consisting  of  only  black  and  white  pixels,  these  colors  have  gray  level  values  of  zero  and  255, 
respectively.  The  MatLab'^  IPT  bwimg  function  generates  an  array  containing  the  binary 
image,  but  the  user  must  define  a  threshold  value  for  it  to  use.  The  threshold  value  is  required 
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by  bwimg  to  create  a  separate  binary  image  from  the  subregion  by  assigning  corresponding 
pixel  gray  level  values  in  the  subregion  -with  greater  than  or  equal  value  to  the  threshold  value 
to  one  in  the  binary  image,  and  those  pixel  values  less  than  the  threshold  value  to  zero  in  the 
binary  image.  For  all  of  the  subregions  analyzed  in  this  thesis,  the  threshold  value  was  set  to 
the  value  of  the  60*  percentile  of  the  gray  level  values  in  the  subregion  under  analysis. 

The  MatLab^  IPT  perimeter  function  operates  on  the  binary  image.  The  perimeter 
function  examines  every  pixel  in  the  binary  image  and  changes  its  value  to  a  one  or  zero  based 
on  the  following  criteria.  For  every  pixel  in  the  binary  image,  the  software  examines  its  four 
nearest  neighbors  and  when  a  pixel  had  at  least  two  neighbors  with  the  values  of  one  and  zero, 
the  software  function  considered  them  to  be  part  of  the  perimeter  of  the  object. 

This  operation  was  effective,  but  as  the  bars  became  smaller  and  the  bar  group  gray 
level  values  became  closer  to  the  background  values,  it  was  necessary  to  add  one  additional 
function  to  eliminate  false  detection  of  perimeter  pixels.  Another  MatLab^  IPT  function  was 
used.  The  majority  filter  function,  called  majority,  examined  every  pixel  in  the  binary  image, 
similarly  to  the  perimeter  function.  However,  instead  of  looking  for  pairs  of  adjacent  ones 
and  zeros,  it  set  each  pixel  value  to  the  value  of  the  most  frequently  occurring  pixel  among  its 
four  neighbors.  This  process  is  analogous  to  giving  each  pixel  a  vote,  and  the  pixels  whose 
vote  represents  the  majority  of  the  vote  determines  the  value  of  the  pixel  under  investigation. 
The  majority  filter  was  applied  twice  to  the  binary  image. 

The  resulting  binary  image  contains  a  mapping  of  the  perimeter  location  from  the  bar 
group  subregion  and  can  be  used  as  a  template  for  use  in  the  construction  of  a  mask  to 
separate  the  bar  group  from  it’s  background.  The  perimeter  mask  produces  a  result  as  shown 
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in  Figure  24.  Figure  24  also  contains  a  plot  of  the  bar  group  cross  section  which  will  be 
discussed  later  in  this  section. 


Figure  24,  Mask  and  Identified  Bar  Group  Cross  Section 


For  smaller  bars,  the  perimeter  mask  does  not  necessarily  isolate  all  three  light  bars  in 
a  bar  group,  but  the  defined  perimeter  is  still  effective  because  the  location  of  the  bar  group  is 
still  demarcated.  This  perimeter  mask  cannot  be  used  to  directly  isolate  the  bar  group  pixels 
from  the  original  image,  however  it  can  be  use  to  construct  the  actual  mask. 

The  actual  mask  is  constructed  directly  from  the  binary  image  of  the  bar  group 
perimeter.  For  this  mask,  the  region  inside  the  perimeter  will  have  pixel  values  set  to  one  and 
the  region  outside  the  perimeter  contains  pixel  values  set  to  zero.  This  type  of  mask  is  often 
called  a  top-hat  filter,  due  to  its  value  of  zero  outside  the  mask  (or  brim)  and  value  of  one 
inside  the  mask  (as  in  the  crown  of  the  top-hat).  This  mask  is  to  be  used  to  isolate  the  bar 
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group  from  the  rest  of  the  selected  region  like  a  cookie  cutter.  It  is  the  same  size  as  the 
original  bar  group  image,  and  the  interior  region  corresponds  in  size  and  location  to  the  region 
in  the  original  image  containing  the  bar  group  to  be  isolated. 

However,  a  few  pieces  of  information  must  be  determined  before  the  mask  can  be 
constructed.  For  instance,  which  pixels  contain  the  top-hat  crown  and  which  contain  the 
brim?  Figure  25  demonstrates  the  method  used  to  determine  the  inner  and  outer  extent  of  the 
bar  group  isolating  mask.  This  figure  depicts  the  binary  image  containing  the  information  on 
the  location  of  the  bar  group  perimeter,  which  is  the  top-hat  crown  location,  within  the 
selected  subimage.  A  horizontal  vector  containing  the  nonzero  pixel  counts  down  each 
column  of  the  binary  image  is  created.  This  process  is  followed  by  the  formation  of  a  vertical 
vector  containing  all  the  nonzero  pixel  counts  across  each  row  of  the  binary  image.  These 
vectors  mark  the  x-range  and  y-range  extents  of  the  binary  mask  perimeter.  A  search  routine 
is  used  to  find  the  x-range  extent  esid  y-range  extent.  The  vector  element  locations  are 
numbered  from  one  to  N,  where  N  is  the  length  of  the  vector.  Within  each  row  and  column 
vector,  it  finds  the  vector  element  locations  containing  nonzero  values’"^.  The  highest  and 
lowest  vector  location  values  are  determined  with  minimum  value  and  maximum  value 
searches,  and  these  addresses  are  saved  as  the  x-range  extent  and  y-range  extent  of  the  mask 
perimeter. 


Due  to  noise,  this  criterion  finds  the  vector  locations  with  fewer  than  two  pixel  values. 
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Figure  25,  Demonstration  of  x-range  and y-range  extent  determination 


The  original  image  is  multiplied  by  this  binary  image,  pixel  by  pixel,  effectively 
removing  the  background  from  the  original  image.  The  resulting  isolated  bar  group  image  is 
placed  into  a  smaller  matrix. 

The  next  two  steps  of  the  block  diagram.  Detect  Bar  Group  Orientation  and  Estimate 
Cross-Section  Amplitudes,  represent  the  crux  of  the  L2  computation.  The  bar  group 
orientation  is  determined  by  using  the  new  isolated  bar  group  image  to  form  horizontal  and 
vertical  vectors  of  pixel  values,  as  was  done  for  the  binary  mask.  Instead  of  containing  the 
sums  of  pixel  values,  these  new  vectors  contain  the  60*  percentile^  ^  values  of  the  respective 
rows  and  columns.  Figure  26  demonstrates  the  formation  of  these  vectors. 


**  The  60*  percentile  is  computed  b>'  finding  the  value  within  a  set  which  is  larger  than  60  percent  of  the  all 
of  the  values  in  the  set;  the  50*  percentile  is  the  median. 
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Figure  26,  60th  Percentile  Raw  and  Column  Vector  Formation 


After  the  application  of  an  exponential  smoothing  filter  [Makxidakis,  Wheelwright,  and 
McGee,  1983:86],  a  second  order  finite  difference  operation  is  performed  on  both  vectors  to 
remove  constant  bias  or  linear  trends.  The  variance  of  both  vectors  is  then  taken.  The  cross 
sections  of  the  bar  groups  are  actually  composed  of  three  adjacent  Bessel  functions  [Gaskill, 
1978:  72].  Since  the  second  derivative,  or  Laplacian,  of  a  Bessel  function  is  another  type  of 
Bessel  function,  the  vector  with  the  largest  variance  is  the  lateral  cross-section  vector  and  the 
orientation  is  defined  accordingly.  Now,  a  mean  lateral  cross-section  vector  can  be  generated 
to  form  an  estimate  of  the  actual  bar  group  cross  section.  Figure  27  shows  the  resulting  mean 
lateral  cross  section  of  the  bar  group  #1  image  selected  earlier.  It  also  displays  the  tenth  order 
polynomial  fit  [MatLab®  Reference  Manual]  created  to  identify  the  pixel  values  and  location 


of  the  maxima  and  minima 


P1F16.TIF,Group#1  .Across  LOF 


Figure  27,  Mean  Lateral  Cross  Section  of  Bar  Group  #I  and  10  th  Order  Polynomial  Fit 


The  L2  value  is  found  by  finding  distances  between  each  peak  and  valley  location  from 
the  mean  lateral  cross  section.  Critical  steps  are  described  in  the  next  two  blocks  of  the  flow 
diagram:  Estimate  Peak  and  Valley  Locations  and  Estimate  Bar  Width.  The  mean  lateral 
cross  section  is  used  to  compute  the  L2  value.  To  start  the  search  for  the  maxima  and  minima 
locations,  the  centroid,  x*,  of  the  polynomial  fit  data  is  computed  according  to  Equation  (2.8). 
This  centroid  value,  and  two  adjacent  data  points,  are  seen  in  Figure  27 


(2.8) 
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This  three  point  comparison  results  in  the  starting  point  for  the  search.  The  starting 
point  for  the  search  is  the  vector  element  containing  the  maximum  value  in  the  polynomial  fit 
data.  The  maximum  value  of  three  points  are  used  to  determine  the  centroid  to  compensate 
for  centroid  location  error  in  the  case  of  poor  selection  of  the  bar  group  subimage.  The 
search  algorithm  is  a  gradient  search  technique  and  is  specifically  written  to  find  the  locations 
of  three  peaks  and  two  valleys  on  the  polynomial  fit  curve.  After  the  search  routine  has  found 
the  three  maxima  and  two  minima,  the  bar  widths  are  computed.  So,  the  smooth  polynomial 
fit  is  used  to  eliminate  false  extrema  in  the  bar  group  cross  section  data  due  to  noise,  but  the 
bar  group  cross  section  data  is  used  to  compute  bar  width  and  modulation.  The  bar  width  is 
estimated  by  finding  the  pixel  differences  in  the  locations  of  adjacent  maxima  and  minima. 

This  is  justified  because  no  edges  exist  to  measure  the  width  of  the  light  and  dark  bars.  The 
original  bars  have  equal  widths.  By  finding  center  to  center  differences  between  adjacent  peak 
and  valley  pixel  locations,  the  actual  width  of  the  bar  group  in  pixels  can  be  estimated. 
Averaging  is  used  to  reduce  errors  introduced  by  under  sampling  and  the  jagged  edges  due  to 
the  noise  introduced  by  the  imaging  process.  Equation  (2.9)  describes  the  estimate  of  the 
actual  bar  width  in  pixels,  Wpixei,  where  the  locations  of  the  maxima  in  Figure  27  were  labeled 
as/7/,  p2,  and ps  from  left  to  right  and  the  minima  locations  were  labeled  v/  and  v^. 

In  order  to  compare  the  estimated  width  to  the  bar  target,  the  width  must  be  converted 
into  units  of  ground  length,  or  centimeters-ground  (cm-ground).  To  make  this  conversion. 


This  conversion  from  gray  level  to  log(£)  units  uses  the  calibration  data  provided  by 
the  user  to  ADiM  during  the  setup  stage.  The  next  two  steps  described  in  the  flow  diagram, 
Convert  From  Gray  Level  to  Exposure  and  Calculate  Bar  Group  Modulation,  were  not  used 
or  required  by  the  OSMPF  to  compute  L2.  This  determination  of  individual  bar  group 
modulations  was  an  attempt  to  determine  another  image  quality  metric  or  reference  and  will 
be  discussed  later  in  this  chapter.  The  method  used  to  ensure  a  good  estimate  of  the 
orientation  and  bar  width  also  provided  data  to  estimate  the  actual  peak  and  valley  amplitude 
values  of  each  bar  group.  The  amplitudes  of  the  three  identified  peak  values  were  averaged 
together  to  form  the  estimate  of  the  bar  group  maximum  amplitude.  The  two  identified  local 
minima,  or  valley  amplitude  values,  were  averaged  to  estimate  the  minimum  amplitude  for  the 
bar  group.  The  maximum  and  minimum  gray  level  amplitudes  were  converted  to  log(£^  units. 
K2  was  calculated  according  to  Equation  (2.5)  and  the  result  was  stored  for  each  bar  group. 

Figure  28  shows  the  results  of  computing  K2  for  every  ADtM  resolved  bar  in  a  bar 
target  image,  for  both  orientations.  The  figure  also  plots  the  average  K2  for  each  bar  width. 
This  data  is  from  the  same  image  as  from  the  full  bar  target  shown  earlier. 
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The  legend  in  the  upper  right  comer  of  this  figure  describes  the  markers  denoting  the 
four  types  of  K2  which  were  computed  to  create  this  graph.  Ten  or  more  measures  of  K2  were 
made  and  recorded  for  each  bar  width  and  in  each  bar  group  orientation.  Please  refer  to 
Figure  1,  in  Chapter  1,  for  a  description  of  the  bar  group  orientations.  The  K2  from  the  bar 
group  measuring  the  in-the-line-of-flight  characteristics  is  called  K2-IOF,  and  the  K2  from  the 
bar  group  measuring  the  across-the-line-of-flighi  characteristics  is  called  K2-XOF.  The 
results  of  these  measures  were  averaged  and  two  more  points  were  plotted  at  each  bar  group 
width.  This  chart  was  an  attempt  to  demonstrate  the  mean  and  variance  statistics.  This  graph 
also  demonstrates  the  consistency  of  ADiM. 

3. 1.1.2.  Finding  K2 

The  second  parameter  used  to  compute  H„in  is  the  modulation  of  the  brightness  panel 
exposure.  The  two  blocks  in  the  flow  diagram.  Select  Region  Containing  Brightness  Panels 
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and  Median  Filter  Brightness  Panels  Suhimage,  begin  the  K2  measurement.  ADiM  parallels 
the  OST  method,  using  the  average  gray  levels  of  the  brightness  panels.  This  value  of  K2  is 
used  to  find  the  ADiM  value  of  Hmin-  The  software  assumes  the  user  can  select  a  sample  of 
the  bright  and  dark  regions  of  the  brightness  panels  within  the  same  rectangular  region.  The 
following  graphic.  Figure  29,  portrays  a  potential  sample  region  for  K2  analysis. 
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Figure  29,  Brightness  Pattel  Selection  For  K2  Analysis 
ADiM  only  differs  from  the  OSMPF  by  the  method  the  dark  and  light  sample  regions 
are  separated.  The  OSMPF  separates  the  two  regions  by  selecting  two  separate  samples,  one 
from  the  dark  brightness  panel  and  one  from  the  light.  The  ADiM  technique  uses  a  single 
subimage  containing  both  regions. 

The  user  starts  ADiM  from  the  initial  screen  shown  in  Figure  22,  and  selects  a  region 
from  the  bar  target  containing  the  brightness  panels.  The  rectangular  area  extracted  must 
contain  pixels  from  both  the  light  and  dark  brightness  panels.  The  selection  is  median  filtered 
to  filter  extraneous  data,  because  it  is  known  that  the  brightness  panels  are  supposed  to 
approximate  a  Lambertian  surface.  Now,  the  subimage  preparation  moves  to  the  next  steps  of 
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the  flow  diagram.  Convert  from  Gray  Level  to  Exposure  and  Calculate  the  Modulation  of  the 
Brightness  Panels. 

An  estimate  of  gray  level  values  for  the  light  and  dark  brightness  panels  is  needed  to 
compute  the  image  modulation,  K^.  Inconsistencies  in  the  images  of  the  brightness  panels 
gave  the  author  cause  to  seek  alternative  ways  to  find  an  estimate  of  K2.  A  histogram  function 
was  used  to  create  a  vector  containing  the  number  of  occurrences,  also  called  counts,  of  the 
gray  levels,  An  example  of  this  histogram  is  shown  in  the  following  figure.  The  consistent 
bimodal  nature  of  the  brightness  panel  data  makes  this  type  of  measurement  possible.  Noise 
effects  are  present  and  can  be  observed  in  the  wide  variances  around  each  of  the  distinct  mean 
gray  level  values. 
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Brightness  Panels  Histogram 


The  number  of  bins  used  to  separate  the  data  was  dependent  on  the  number  of 
samples.  In  general,  for  large  sample  sizes,  the  number  of  bins  in  the  histogram  were 
computed  as  =  A:  log(n),  where  is  the  number  bins,  k  is  a  constant  with  values  between  2 
and  5,  and  n  is  the  number  of  samples  available.  The  MatLab®  Histogram  function,  hist, 
returns  two  vectors,  a  vector  containing  the  gray  level  counts  for  each  bin  and  a  vector 
containing  the  center  value  of  each  bin.  The  bins  containing  the  greatest  frequencies  of  gray 
level  values  are  found  for  both  the  high  (light)  and  the  low  (dark)  gray  level  values.  The 
values  are  then  averaged,  converted  to  log(£)  units,  and  used  to  calculate  K2. 
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3.2.  Discussion  of  Data  Specifics 


The  idiosyncrasies  associated  with  the  data  drove  the  development  of  ADiM  from  the 
beginning.  The  high  noise  level  and  the  undersampling  present  made  creating  an  automatic 
solution  to  the  processing  needs  of  the  OSMPF  very  difficult.  Some  of  the  adaptations  made 
to  ADiM  to  handle  these  properties  have  been  described  earlier,  but  in  the  next  sections, 
additional  topics  will  be  covered. 

Each  system  in  the  process  to  satisfy  the  OST  calibration  requirements  introduced 
uncertainty.  Film  is  a  nonlinear  detector  and  adds  uncertainty  to  the  process  if  it  is  not  used 
properly.  The  digitization  of  the  film  data  introduces  other  problems  with  sampling  size.  The 
effect  of  the  imaging  process  on  bar  target  edges  has  been  shown,  but  additional  information 
on  finding  thresholds  to  demarcate  these  objects  is  included  below. 

3.2.1.  Noise  Effects  in  the  Data 

In  truth,  a  bar  group  is  really  not  composed  of  light  and  dark  bars.  Each  bar  group 
effectively  consists  of  white  bars  painted  on  a  black  background  and  so  the  set  of  contours 
around  the  artifacts  with  the  strongest  signal  return  from  the  ground  target,  and  if  the 
selection  of  the  threshold  value  is  suitable,  the  perimeter  gray  level  will  be  larger  than  all  of 
the  background  gray  levels,  such  that  only  the  bar  group  signal  has  gray  level  values  which 
exceed  the  perimeter  value. 

To  understand  the  importance  and  difficulty  of  finding  a  thresholding  value,  the 
stochastic  nature  of  the  data  must  be  discussed.  Although  all  of  the  images  considered  in  this 
thesis  have  pixel  values  between  1  and  256,  the  variance  and  mean  of  the  data  is  not 
predictable  between  different  images,  nor  are  the  statistics  of  the  gray  level  values  consistent 
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within  a  given  calibration  target  image  from  bar  group  to  bar  group.  Few  assumptions  can  be 
made  about  the  distributions  of  the  gray  level  values  in  any  sample  image,  because  the 
probabilistic  distributions  of  gray  levels  in  all  images  are  asymmetric,  bimodal,  and  therefore 
non-Gaussian  [Kreyszig,  1983;  981], 

At  the  end  of  Chapter  2  a  discussion  of  the  lack  of  homogeneity  in  the  brightness 
panels  demonstrated  the  inconsistency  and  roughness  of  the  gray  level  values  across  larger 
surfaces  on  each  target  with  apparently  the  same  color  of  reflective  paint.  The  variance  in 
their  surface  gray  level  values  appears  to  be  normally  distributed,  but  is  very  wide.  If  the  bar 
group  sampled  from  the  target  image  contains  too  many  samples  from  the  surrounding  area, 
the  inconsistencies  bias  the  ADiM  calculations  and  the  perimeter  edges  of  the  light  and  dark 
bars  cannot  be  distinguished. 

Figure  3 1  portrays  a  histogram  of  the  gray  levels  from  an  image  of  the  WPAFB  target 
taken  on  25  Aug  95.  Notice  the  two  distinct  modes  with  local  minima  at  the  left,  right,  and 
center  of  the  distribution.  In  many  image  processing  applications,  edges  can  be  identified  by 
choosing  the  threshold  value  based  on  the  histogram  minimum  value  [Castleman,  1996:247]. 
Using  the  histogram  central  minima  present  in  the  histogram,  one  could  infer  that  the  bar 
group  edges  are  defined  by  the  transition  from  dark  to  light  values  around  the  pixel  gray  level 
of  160.  So,  if  each  bar  group  had  enough  samples  in  respective  subimages  to  form  a 
meaningful  histogram,  and  all  of  the  bar  groups  appeared  to  have  bimodal  gray  level 
distributions  with  the  central  minima  occurring  around  the  160  gray  level  range,  determining 
the  edge  locations  among  each  of  the  bar  triads  would  be  trivial.  The  histogram  of  such  a 
target  would  show  how  to  threshold  gray  level  values  in  order  to  differentiate  signal  from 
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noise;  unfortunately,  it  does  not  work  for  the  situation  defined  by  the  OS  Treaty.  The 
proportionally  shrinking  bar  group  dimensions  on  the  OS  target  have  correspondingly 
decreasing  illuminance  amplitudes.  This  decreasing  amplitude  is  a  decrease  in  signal  level  and 
it  is  the  result  of  the  convolution  of  the  optical  system  MTF  and  the  reflected  ground  image  of 
the  target.  The  physics  of  the  effects  of  the  optical  system  MTF  were  discussed  in  more  detail 
in  Chapter  2. 

For  now,  the  important  fact  is  that  the  gray  level  values  along  the  perimeters  of  the 
light  and  dark  bars  on  the  OS  target  decrease  as  the  bar  group  width  and  length  dimensions 
decrease,  and  the  number  of  samples  within  each  subimage  decreases  as  each  group  gets 
smaller  in  order  to  reduce  the  introduction  of  noise  from  artifacts  other  than  bar  groups  which 
are  present  in  every  OS  target.  Notice  the  60***  percentile  of  this  data  is  near  the  gray  level 
value  of  160. 
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P1F16  Target  Data 


But  in  the  case  of  the  OS  data,  the  bimodal  nature  of  the  data  cannot  be  used  because 
the  distribution  of  gray  levels  changes  and  not  enough  samples  exist  in  the  smaller  bar  groups, 
as  shown  in  Figure  32,  Almost  all  of  the  analyzed  target  data  sets  in  this  thesis  had  GRD 
values  less  that  30  cm,  which  is  the  bar  width  of  bar  group  #9.  The  details  available  in  a 
histogram  depend  on  the  number  of  samples  used  in  its  computation;  since  the  following 
figure  represents  the  histogram  computed  from  a  resolved  bar  group,  its  lack  of  available 
structure  proves  the  inappropriateness  for  use  in  finding  a  proper  threshold  value. 
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BarGroLp#9,P1F16 


Figure  32,  Histogram  of  Bar  Group  #9 


3.2.2.  For  smaller  bars,  poor  choice  in  bar  group  location  makes  automatic 
peak  and  valley  identification  difficult 

As  the  bar  widths  narrow,  the  amplitudes  of  the  peaks  and  valleys  decrease,  and  the 
relative  distances  between  the  peaks  and  valleys  also  decrease.  This  is  a  result  of  the  system 
MTF  as  discussed  in  Chapter  2.  In  other  words,  as  the  width  of  the  bar  groups  get  smaller, 
the  peaks  and  valleys  become  more  and  more  difficult  to  distinguish,  until  only  a  single  peak 
remains.  The  smallest  bar  group  for  which  the  peaks  and  valleys  are  discernible  should  be  the 
GRD,  but  it  is  possible  for  noise  or  atmospheric  distortion  to  blur  a  bar  group  image  beyond 
the  true  system  GRD.  The  resulting  GRD  for  the  target  could  be  greater  than  or  less  than  the 
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noise  free  result.  Bar  group  #9  is  significantly  smaller  than  bar  group  #1,  and  is  shown  below 
in  Figure  33. 


Figure  33,  Example  of  the  Effect  of  Decreasing  Barwidth  on  Peak  and  Valley  Amplitudes  of 
Bar  Group  #9  (30  cm-ground  bar  width,  no  filtering) 

The  edges  of  the  bar  groups  are  rounded  through  the  imaging  process.  For  all  bar 
groups,  MTF  and  noise  effects  created  jagged  regions  along  each  bar  length  in  the  group,  but 
in  a  small  bright  bar  these  jagged  regions  represented  amplitude  changes  which  could  be  as 
large  as  the  amplitude  of  the  bar  group  itself  The  60*  percentile  threshold  value  allowed 
ADiM  to  define  a  legitimate  contour  around  a  group  of  interest  and  the  MatLab®  IPT 
functions  provided  a  means  to  extract  the  most  likely  contours  inside  of  the  noisy  data.  For 
small  bars,  it  was  still  often  necessary  to  try  several  different  subimage  selections  before 
ADiM  successfully  identified  the  bar  group. 


79 


If  the  OSMPF  had  to  calculate  the  K2  value  from  the  bar  group  bounding  the  system 
GRD  and  not  from  the  brightness  panels,  the  issue  of  bar  group  edge  selection  would  become 
a  significant  issue,  because  the  K2  value  calculated  from  the  bar  group  varied  up  to  25% 
depending  on  the  ROI  selection.  The  K2  estimate  is  reduced  by  the  imaging  process  edge 
smoothing,  because  it  reduces  the  average  maximum  amplitude  measured.  Then  the  mean 
lateral  cross  section  is  computed,  and  the  amplitudes  of  the  maxima  are  underestimated. 

ADiM  avoids  this  by  cropping  the  pixels  from  the  edges  before  computing  the  mean  lateral 
cross  section.  The  figure  below  shows  the  reduced  peak  amplitude  at  the  narrow  edge  of  the 
bar  group  image  By  increasing  the  sample  rate,  without  increasing  the  noise  of  the  system, 
this  sensitivity  could  be  reduced  by  increasing  the  knowledge  of  the  region  around  the  smallest 
bar  groups. 
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Figure  34,  Gray  Level  hjtensity  Image  of  Bar  Group  #1(75. 6  cm-ground  har  width) 
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3.2.1.  Photographic  Film  Usage 

The  use  of  photographic  film  is  well  documented  [Slater,  1983],  Photographic  film 
consists  of  two  parts,  the  photosensitive  emulsion  which  reacts  to  exposure  to  light  and  the 
acetate  to  which  it  is  adhered.  For  hlack  and  white  photography,  this  emulsion  is  a  silver 
nitrate  based  compound.  When  the  emulsion  is  exposed  to  light  and  put  through  the  liquid 
chemical  developing  process,  the  parts  of  the  emulsion  exposed  to  light  remain  on  the  acetate 
in  proportion  to  the  amount  of  light  to  which  it  was  exposed.  At  the  end  of  the  developing 
process  the  emulsion  remaining  on  the  film  has  an  optical  density  which  is  linearly  related  to 
the  logarithm  of  the  exposure  from  the  incident  light.  The  image  present  in  the  incident  light 
has  its  darker  details  recorded  in  thinner  regions  of  emulsion  and  its  brighter  regions  recorded 
in  thicker  regions.  Lower  values  of  optical  density  (OD)  correspond  to  a  thinner  presence  of 
emulsion,  and  higher  exposure  level  to  light.  The  plot  of  OD  versus  log£  is  called  the 
Hurter-Driffield  (H-D)  curve.  Although  not  shown  here,  the  curve  is  an  increasing  s-shaped 
curve  and  a  similar  curve  is  discussed  below. 

For  the  OST,  if  the  film  is  exposed  to  too  much  light  or  not  enough  light,  the 
information  contained  on  the  film  will  be  seriously  degraded  or  misinterpreted.  The  OSMPF 
have  two  sources  of  information  to  ensure  the  film  is  used  in  the  proper  range  of  exposure. 
NBS  film  strips  of  known  optical  density  allow  the  OSMPF  to  calibrate  their  optical  scanners 
from  gray  level  values  to  optical  density.  They  also  have  an  instrument  used  on  each  roll  of 
film  used  for  an  OS  Treaty  flight  or  practice  flight  which  exposes  the  mission  film  beginning 
(head)  and  end  (tail)  to  a  series  of  know  exposures,  in  order  to  convert  the  digital  scanner 
output  from  gray  levels  to  units  of  the  logarithm  of  exposure,  \o%E  .  The  logJ?  values  from 
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the  brightness  panels  are  computed  and  plotted  on  the  H-D  curve  for  the  film.  As  long  as  the 
plotted  values  reside  in  the  linear  region  of  this  curve,  the  OSMPF  considers  the  film  properly 
exposed. 

The  data  analyzed  in  this  thesis  had  been  analyzed  by  the  OSMPF  and  was  known  to 
originate  fi'om  properly  exposed  film,  so  there  was  no  need  to  generate  a  calibration  curve  of 
optical  density  of  the  film  at  different  light  levels  to  the  H-D  curve.  The  conversion  for  the 
KODAK  5057  scanner  output  from  gray  level  to  logF^  was  made  in  ADiM,  The  shape  of 
this  curve  is  the  same  shape  as  the  H-D  curve,  because  the  scanner  gray  level  output  is  linearly 
proportional  to  the  OD  of  the  input  photograph.  A  minimum  sum  of  squares  criterion  is  used 
to  determine  the  data  points  to  be  used  in  the  computation  of  the  linear  model  to  describe  the 
gray  level  and  log(JE')  relationship.  This  model  is  used  to  predict  the  logCE')  values  from  the 
gray  level  values  and  fits  to  the  central  linear  region  of  the  elongated  s-shaped  curve  [Neter, 
Wasserman,  and  Kunter,  1990:26]. 

For  this  criterion,  the  centroid  of  the  region  is  found  and  pairs  data  points  smaller  and 
larger  than  the  centroid  are  added  to  the  model  iteratively.  With  each  pair  of  additional  data 
points,  the  regression  coefficients  were  computed,  the  mean  square  error  (MSB)  was 
computed  between  the  data  and  the  model  prediction.  To  select  the  final  model,  the  MSB 
values  were  compared  for  each.  The  model  containing  more  than  three  points  and  the 
smallest  MSB  was  used  for  the  calibration  curve.  Figure  35  shows  the  linear  region  and  the 
distinctive  s-shape.  A  plot  of  the  log(£')  versus  the  gray  levels  generated  by  the  KODAK 
5057  image  scanner  are  shown  below. 
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Figure  35,  Log  E  V5.  2 1 -Step  Wedge  Gray  Level  Curve  Demonstration  of  the  S-Shape  and 

Linear  Region 

K2  is  determined  with  an  estimate  of  the  gray  level  values  from  the  light  and  dark 
brightness  panels,  which  are  converted  into  values  with  units  of  exposure  (Lumens  * 
integration  time).  The  logE'  values  for  the  three  maxima  are  averaged  and  the  values  for  the 
two  minima  are  averaged.  The  difference  is  taken  and  used  in  Equation  (2.4)  and  then  in 
Equation  (2.5)  in  order  to  compute  the  actual  value  for  K2  to  be  used  in  the  calculation  for 

H rnin . 

The  complete  data  sets  in  this  thesis  were  analyzed  using  data  digitized  by  the 
KODAK  5057  scanner.  One  incomplete  set  of  data  from  the  Perkin-Elmer  Digital  Scanner 
was  considered,  and  the  ADiM  software  was  able  to  resolve  bar  groups  commensurate  with 
OSMPF  results. 
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3.3.  ADiM  Development  Further  Elaboration 


ADiM  does  not  re-implement  the  OS  method  for  finding  the  GRD.  In  fact,  it  cannot 
due  to  the  excessive  variance  and  noise  in  the  data,  as  well  as,  the  undersampling  of  the 
smaller  bar  groups  of  interest.  ADiM  increases  the  confidence  in  estimates  of  the  bar  group 
signal  amplitudes  and  locations  in  order  to  improve  the  estimates  of  L2  and  K2.  To  do  this, 
ADiM,  like  most  digital  techniques,  assumes  symmetric  error  distributions  between  the  pixel 
values  and  locations  of  the  recorded  bar  group  and  those  corresponding  to  them  in  the 
theoretical  bar  group  image.  As  a  result  of  this  assumption,  it  can  take  advantage  of  the 
variance  reduction  achieved  through  averaging  related  values  in  the  raw  data.  Averaging 
reduces  the  effects  of  symmetrically  distributed,  random,  unbiased,  noise. 

In  general,  averaging  is  applied  within  a  single  data  set  or  image  or  with  an  ensemble 
of  repeated  measurements  of  a  data  source  such  as  series  of  images.  The  type  of  data 
currently  taken  for  the  OS  regime  does  not  lend  itself  to  ensemble  averaging  techniques  for 
several  reasons.  Since  the  frame  rate  is  only  fast  enough  to  minimize  blurring  due  to  aircraft 
motion  and  the  aircraft  flies  at  different  altitudes  and  attitudes  for  each  pass  over  the  target. 

In  addition,  the  number  of  flights  required  for  each  configuration  to  generate  a  significant 
ensemble  is  too  expensive.  Hence,  ADiM  is  applied  to  individual  target  subregions  and  uses 
row  and  column  averaging  to  estimate  the  bar  group  width  and  amplitude. 

Median  filters  and  lowpass  filters  are  widely  used  to  reduce  and  smooth  random  noise 
effects  in  data  and  thereby  improve  the  estimated  values  of  signal  amplitudes  [Castleman, 
1996:  207,  247].  Both  filters  are  used  by  ADiM  to  calculate  the  two  components  of 
The  first  operation  performed  toward  the  determination  of  the  GRD  on  the  selected  subregion 
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of  the  image  is  a  median  filter.  Each  pixel  is  replaced  with  the  median,  or  50*  percentile, 
value  among  it  and  its  eight  adjacent  neighbors.  The  median  filter  displaces  “salt  and  pepper” 
type  noise  pixels  toward  the  edges  of  the  image.  This  noise  is  caused  by  pixels  with  gray  level 
values  significantly  different  from  their  surrounding  neighbors.  Since  the  calibration  target  is 
known  not  to  have  any  intentional  high  or  low  reflective  speckle  associated  with  its  surface, 
this  median  filter  is  not  removing  any  usefiil  information.  It  tends  to  sharpen  edges,  but  is 
non-linear  and  non-reversible  [Castleman,  1996:247,  Thompson,  1995:79,83]. 

After  the  first  binary  image  is  formed,  the  produces  a  binary  edge  map  image  of  the 
locations  of  the  bar  group  perimeter  pixels.  The  edges  are  assigned  the  value  of  one,  and  the 
remaining  pixels  are  set  to  values  of  zero.  The  edges  that  are  found  are  defined  corresponding 
to  the  edge  location  in  the  actual  image.  IheMatLab®  IPT  is  used  to  eliminate  the  pixels 
which  were  spuriously  assigned  the  value  one.  Each  pixel  is  replaced  by  the  value  of  the 
majority  value  of  its  four  adjacent  pixels,  not  including  diagonally  adjacent  pixel  location.  Due 
to  these  MatLab®  IPT  functions  and  the  median  filter,  very  few  single  pixels  will  have  the 
value  of  one  outside  the  bar  group  perimeter  in  the  binary  edge  map  image. 
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rv.  Results 


In  this  chapter,  the  results  of  the  data  processing  are  analyzed.  With  only  one 
exception,  the  data  comes  from  Open  Skies  Certification  Mission  95-MOC-007,  25-27 
August  1995.  One  set  from  another  mission,  95-MOC-006  on  25  April  1995  is  used.  The 
former  mission  used  the  WPAFB  target  for  calibration  and  the  latter  used  the  51-51  Target  at 
Eglin  AFB,  Florida.  The  data  is  catalogued  first  by  the  mission  name  and  the  sensor  used  to 
collect  it,  then  by  the  film  type  used,  followed  by  any  filters  used,  and  finally  by  the  physical 
pass  number  over  the  target  and  the  frame  of  film  on  which  the  data  is  recorded.  All  data  was 
taken  with  the  OS  KS-87  sensor,  which  is  a  framing  camera  with  an  fr4.5,  76.2  mm  focal 
length  lens.  The  KS-87  is  a  single  camera  system  and  can  be  mounted  in  three  separate 
positions  on  the  OS  OC-135  aircraft.  With  these  positions  it  can  record  data  looking 
vertically  down  to  the  ground,  or  out  two  different  aircraft  fuselage  portals  to  cover  the  land 
areas  along  the  left  and  right  sides  of  the  aircraft.  The  data  analysis  in  this  thesis  only 
considers  the  data  from  the  vertical  and  left  looking  positions.  The  MTF  measurements  of  the 
KS-87  lens  system  are  included  in  Appendix  2  for  several  incident  angles  and  positions.  On- 
axis  performance  is  near  the  dififaction  limit,  but  as  the  incident  angle  diverges  from  the  center 
line,  the  performance  degrades. 

The  data  in  this  thesis  was  collected  at  altitudes  which  vary  from  1460  m  to  7010  m 
above  ground  level.  The  weather  conditions  were  not  considered  to  have  a  strong  influence 
on  the  collection.  Each  configuration  for  sensor,  position,  filter,  and  film  must  be  flown  over 
the  target  of  interest,  as  many  times  as  is  agreed  to  between  the  aerial  imagery  collectors  and 
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the  land  area  governors.  In  this  thesis,  five  passes  in  four  configurations  are  explored.  An 
additional  set  of  data  fi'om  a  separate  mission  over  the  Eglin  AFB  5151  is  included  as  an 
example  of  the  effectiveness  of  the  ADiM  to  find  L2  for  added  scenarios.  The  single  set  of 
EAFB  target  data  analyzed  is  exhibited  without  any  information  about  the  modulation  or  film 
characteristics  needed  to  calculate  Hmm,  but  it  is  valuable  in  that  it  demonstrates  the  ADiM  is 
stable  for  cases  from  different  environments.  Some  analysis  of  simulated  data  is  performed — 
not  to  calculate  Hmm,  but  to  demonstrate  that  the  success  of  the  order  statistics  method  used 
on  actual  target  images  is  not  by  chance,  but  is  an  analytically  sound  method  for  estimating  the 
values  for  the  bar  group  image  extrema.  Each  image  from  a  physical  pass  over  a  target  was 
recorded  on  a  frame  of  film  and  then  the  target  region  was  digitized  as  described  previously 
and  converted  directly  to  a  TIFF  file.  The  scanner  employed  here  is  capable  of  producing 
images  of  several  Megabytes  in  volume,  so  it  is  usually  necessary  to  crop  all  of  the 
surrounding  area  around  the  target  and  reduce  the  file  size  to  hundreds  of  kilobytes. 

In  an  ideal  experiment,  this  TIFF  file  would  contain  a  header  containing  useful 
information  about  the  image  and  the  binary  representation  of  each  sample  pixel  from  the 
image.  No  processing  or  averaging  would  be  done  at  this  stage,  but  much  of  the  unnecessary 
imagery  information  would  be  cropped  from  the  area  surrounding  the  target.  Unfortunately, 
one  of  the  four  sets  of  data  received  from  the  OSMPF  had  accidentally  been  digitally 
resampled  while  using  Adobe®  Photoshop®  software  for  cropping  the  data.  The  sampling  size 
was  decreased  (the  sampling  rate  was  increased)  within  Photoshop®  by  interpolating  between 
the  raw  image  data  pixels  and  adding  additional  points  between  each  pair.  This  causes  an 
averaging  effect  which  is  similar  to  low  pass  filtering.  One  of  these  interpolated  sets  of  target 
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data  was  a  duplication  of  a  previously  provided  data  set.  It  will  be  used  to  show  the  effect  of 
the  resampling  of  the  data  is  minimal  and  will  not  prevent  the  use  of  this  data  for  analysis  in 
this  thesis.  This  result  should  not  be  surprising  because  the  nature  of  the  method  described 
here  depends  heavily  on  averaging. 

Using  the  techniques  of  ADiM,  the  results  that  follow  were  found.  In  the  first  set  of 
results,  the  OSMPF  values  for  L2,  K2,  and  Hmm  are  compared  directly  to  the  ADiM  values. 
Due  to  a  consistent  residual  between  the  OSMPF  value  for  K2  and  the  value  for  K2  found  in 
this  thesis,  the  next  results  use  the  OSMPF  values  for  K2  to  calculate  Hmin  for  the  L2  values 
from  both  techniques.  Although  more  research  and  testing  needs  to  be  performed,  this 
difference  in  K2  values  is  attributed  to  the  lack  of  NBS  calibration  for  the  KODAK  5057 
scanner.  The  Zeiss  microdensitometer  is  regularly  calibrated  in  accordance  with  NBS 
guidelines. 

4.1.  Result  Comparison 

Figure  36  shows  a  scatter  plot  of//m,„  values  determined  by  the  OSMPF  (horizontal 
axis)  and  the  Hmm  values  determined  by  ADiM  (vertical  axis).  It  supports  the  supposition  that 
the  results  from  both  techniques  are  related  in  a  simple  linear  fashion.  The  result  of  this 
supposition  is  that  each  calculation  can  be  predicted  by  knowing  the  other.  Finding  the  linear 
fit  for  the  two  sets  of  data  containing  25  data  points  in  the  form  of  a  simple  linear  model  yields 
a  correlation  coefficient  of  0.886.  This  means  almost  89  percent  of  the  variation  in  the 
OSMPF  Hmin  value  is  explained  by  the  variation  in  the  ADiM  value.  The  F-Statistic  for  the 
model  was  1.96E-08.  The  p-valnes  for  the  estimated  slope  and  intercept  coefficients,  pOmA 
pi,  were  0.21872  and  1.96E-08,  respectively.  Thus,  the  parameters  and  model  are 
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statistically  significant.  This  demonstrates  that  ADiM  values  are  correlated  to  OSMPF,  and 
the  methods  changed  in  order  to  adapt  the  OSMPF  techniques  to  digital  processing  are 
consistent  with  the  original  method. 


Hmin  for  OSMPF  vs  Hmin  for  ADiM 


OSMPF  Hmin  (m) 

Figure  36,  Graph  demonstrating  a  linear  correlation  between  OSMPF  and  ADiM 

measurements 


Explained  variation  does  add  to  the  credibility  of  the  results,  but  a  scatter  plot  of  both 
OSMPF  and  ADiM  computed  results  of  Hmi„  in  Figures  36  and  37  show  that  none  of  the 
ADiM  results  exactly  matched  the  OSMPF  values  for  H„i„ ,  so  some  systematic  differences  in 
the  two  methods  exist.  This  difference  in  results  can  be  traced  most  often  to  different  values 
of  the  calculated  film  exposure  modulation  values  between  the  two  measurements. 
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Hmin  from  OSMPF  and  ADiM 
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Figure  37,  Graph  shows  a  relationship  between  OSMPF  and  ADiM  Hmm  measurements 

Figure  38  displays  the  same  data  as  Figure  37,  but  the  order  of  the  data  points  has 
been  sorted  by  the  increasing  OSMPF  H„i„  values.  The  oscillation  of  the  ADiM  results  for 
data  points  numbered  10-20  could  support  the  presence  of  some  systematic  error  in  ADiM. 
Although,  it  must  be  said  that  an  difference  in  L2  value  of  even  one  bar  could  result  in  such 
discrepancies. 
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Hmin  Values  for  OSMPF  and  ADiM 


Corresponding  Sorted  Data  Points 

Figure  38,  Same  data  as  Fig  33 — except  the  values  were  sorted  in  ascending  order 

A  simple  linear  model  was  built  with  the  least  squares  method  to  determine  the 
correlation  between  the  ADiM  and  OSMPF  results.  In  order  to  assess  the  value  of  a 
simple  linear  model,  a  diagnostic  plot  of  the  residuals,  the  difference  between  the  estimated 
values  of  from  the  model  and  the  actual  Hmin  values  calculated  by  the  OSMPF,  versus  the 
calculated  OSMPF  Hmm  values  are  shown  in  Figure  39  [Neter,  Wasserman,  Kunter  1990: 117], 
The  left  axis  measures  the  aforementioned  residuals  between  the  model  and  the  OSMPF 
values.  The  three  symbol  types  correspond  to  three  different  film  types  used  in  the  analysis. 
This  structure  could  imply  that  higher  order  terms  are  needed  in  the  linear  model  to  explain 
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additional  variance  in  the  OSMPF  results  with  the  ADiM  results.  Note  how  the  types  of  film 
group  the  residuals  into  three  distinct  regions. 


OSMPF  Hmin  Values  (m) 

Figure  39,  Linear  Model  Residuals  (Error)  plotted  against  OSMPF  Hmin 


Another  parameter  which  varied  within  the  data  analyzed  for  this  thesis  was  film  type. 
Figure  39  demonstrates  the  possibility  that  different  films  affect  the  residuals  between 
measured  and  modeled  values.  This  structure  in  the  residual  plot  could  also  imply  that 
improper  calibration  data  was  used  to  compute  K2,  since  a  sample  of  each  type  of  film  must  be 
used  in  order  to  compute  the  modulation  and  translate  the  digitized  values  into  measures  of 
luminance.  Nevertheless,  with  linearity  of  Figure  36  and  the  similar  trends  between  the  ADiM 
and  OSMPF  values  for  all  of  the  data  runs ,  the  simple  linear  model  will  be  accepted  as  valid. 
The  linearity  of  the  comparison  between  the  two  sources  of  computed  Hmm,  albeit  under  poor 
calibration  conditions,  supports  the  validity  of  the  technique. 
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V.  Conclusion 


Statistics  and  commercial  image  processing  tools  can  be  used  to  develop  a  less 
subjective  method  of  determining  the  GRD  and  K2  of  an  OST  ground  target.  The  method 
developed  in  this  thesis  is  not  fully  automated,  and  a  human  is  used  to  select  a  region  in  which 
it  is  knoAvn  a  bar  group  exists.  Undersampling  of  the  bar  groups  made  the  creation  of  a  fully 
automated  technique  too  formidable. 

In  fact,  the  most  promising  option,  in  the  opinion  of  the  author,  to  automate  the 
OSMPF  process  involved  fitting  a  digital  template  over  the  digitized  image  and  scanning  these 
predefined  areas  region  by  region  for  bar  groups  and  brightness  panels  to  compute  the  GRD 
and  K2.  This  template  concept  could  work  with  the  current  digitized  resolution,  but  obtaining 
an  accurate  fit  would  be  extremely  difficult  for  the  undersampled  regions,  whose  widths  are 
often  representative  of  the  GRD.  This  could  be  solved  if  the  calibration  flights  were  flown  at 
lower  altitudes,  or  a  higher  resolution  digitizer  could  be  used.  Thus,  it  is  possible  the 
techniques  developed  in  this  thesis  to  locate  and  identify  the  bar  group  extent  and  orientation 
could  be  extended  to  remove  the  human  in  the  loop,  but  this  step  is  left  for  the  next 
investigator  to  uncover. 

The  sample  size  of  the  OST  target  with  the  photographic  film  would  have  greatly 
simplified  this  task,  because  it  was  on  the  order  of  single  micrometers,  which  translates  to 
around  20  mm  of  ground  dimension  for  the  calibration  runs  flown  under  2000  m. 
Unfortunately,  the  KODAK  5057  digitizer  increased  the  sample  size  three  times  for  the  same 
data  runs.  Due  to  the  large  sample  size  introduced  with  the  digitizer,  the  smaller  bars  within 
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each  target  were  often  only  defined  by  two  or  three  samples  across  their  width,  especially  for 
the  image  data  collected  at  altitudes  over  1800  m.  Undersampled  images  result  in  hard  to 
quantify  GRD,  because  the  edges  of  the  bar  groups  are  difficult  to  isolate.  However,  the 
shape  and  location  of  the  bars  were  known  a  priori  and  this  knowledge  made  identification 
much  easier. 

In  order  to  fijlly  automate  this  technique,  higher  sample  rates  should  be  used.  A  higher 
sampling  rate  would  provide  more  samples  which  could  be  used  to  generate  clearly  bimodal 
histograms.  These  histograms  could  be  interpolated  and  used  to  determine  the  threshold 
values  defining  the  transition  from  light  to  dark  bars.  The  OST  participants  could  experiment 
with  flying  calibration  missions  at  lower  altitudes  to  increase  sample  rate. 

The  specialized  tool  developed  in  this  thesis  is  fianctional,  but  better  results  would  be 
reaped  from  a  change  in  data  collection  technique.  The  current  technique  assumes  the 
participants  will  only  use  single  frames  of  images  of  double  or  triple  bar  targets  for  calibration 
analysis,  but  with  the  cost  of  computing  power  decreasing,  good  confidence  in  the 
identification  of  details  much  smaller  than  those  available  on  a  single  frame  can  be  achieved. 
Faster  frame  rates  allow  ensemble  averaging  of  the  data  and  reduce  the  normally  distributed 
error  introduced  to  the  images  by  atmospheric,  optical,  and  electronic  noise. 

An  even  further  change  in  the  data  collection  paradigm  involves  a  radical  change  in  the 
ground  target  itself— to  a  configuration  allowing  a  continuous  comparison  of  identified  bar 
edge  distinctions  with  decreasing  bar  width  dimension.  A  large  black  and  white  “spoke 
target”  [Goodman,  1988: 126]  could  work  properly  and  the  analysis  would  use  the  diameter 
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cross-sections  for  averaging  to  produce  a  graph  depicting  the  smallest  angular  separation 
achieved  in  a  given  photograph. 

Another  improvement  in  the  OSMPF  method  for  calibration  would  be  to  increase  the 
number  of  samples  in  the  K2  computation.  Currently,  they  use  a  10  pixel  by  10  pixel  region. 
Due  to  the  lack  of  homogeneity  within  the  brightness  panels,  more  samples  are  useful  in 
increasing  the  confidence  in  this  computation.  This  improvement  could  be  enacted  with  very 
little  change  in  current  procedure. 
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Appendix  I  Numerical  Results 


Table^^^NumericalResults^omparin^^^OSA^FAnal^sisJo^DiA^ 


Target 

Sensor 

Film 

Altitude 

(meters) 

K2  (cm) 

OSMPF 

K2  (cm) 

ADiM 

L2  (cm) 

OSMPF 

80%  rule 

L2  (cm) 

ADiM 

lOF/XOF 

Hmin 

OSMPF 

(meters) 

F/inirt 

ADiM 

(meters) 

EAFB 

KS-87 

3404 

7010 

na 

na 

#25/#21 

#25/#21 

na 

na 

WPAFB 

KS-87V 

3404 

1697 

0.627 

0.5159 

23.8 

16.1 

1747.4 

1698.8 

1666 

0.641 

0.5136 

26.7 

23.8/” 

1512.3 

1812.7 

1630 

0.634 

0.5143 

26.7 

23,8/” 

1486.95 

1834.5 

<5 

1610 

0.641 

0.4864 

23.8 

23.8/” 

1640.7 

1857.6 

<  ? 

1608 

0.648 

0.4854 

26.7 

26,7/23.8 

1452.5 

1755.7 

WPAFB 

KS-87V 

3404 

1697 

0.627 

0.5113 

21.2 

23.8/” 

1747.4 

2032.8 

c  > 

(  ) 

1666 

0.641 

0.5135 

21.2 

21.2/” 

1512.3 

2106.9 

C  ) 

«  ) 

1630 

0.634 

0.5054 

23.81 

21.2/” 

1486.95 

2076.2 

c  ? 

4  7 

1610 

0.641 

0.518 

21.2 

21.2/18.9 

1640.7 

2151.5 

c ; 

«  ’ 

1608 

0.648 

0.5052 

23.8 

21.2/18.9 

1452.5 

2173.2 

WPAFB 

KS-87L 

3412 

2746 

0.511 

0.5063 

21.2 

30/” 

3478.7 

2469.7 

i  » 

t  j 

2743 

0.502 

0.4944 

21.2 

30/” 

3502.8 

2493.6 

( ) 

2741 

0.511 

0.489 

23.8 

30/” 

3093.2 

2504.1 

( ) 

2746 

0,519 

0.4839 

21,2 

26.7/30 

3481.2 

2678.2 

«  3 

2741 

0.519 

0.4984 

23.8 

30/” 

3071.6 

2482.7 

WBAFB 

KS-87V 

sao50 

1467 

0.536 

0.5039 

33.7 

23.9 

1145.8 

26.7 

(  ) 

1470 

0.552 

0.4996 

30.0 

23.9 

1271.7 

23.9 

1472 

0.568 

0.5670 

33,7 

23.9 

1120.1 

21.2 

1469 

0.56 

0.5280 

33.7 

26.7 

1125.0 

21.2 

1484 

0.576 

0.5186 

33.7 

26.7 

1107.0 

21.2 

WPAFB 

KS-87L 

SO-050 

1484 

0.654 

0.4003 

33.7 

30 

1059.9 

30 

‘  ’ 

1469 

0,641 

0.3181 

37.8 

26.7 

943.2 

26.7 

1469 

0.646 

0.4056 

33.7 

26.7 

1054.9 

26.7 

'  *’ 

1468 

0.661 

0.4159 

30 

26.7 

1171.0 

26.7 

1470 

0.627 

0.4050 

33,67 

33.6 

1069.9 

33.6 

Appendix  II  MatLab  ADiM  Code 


ADiM  Source  Code.  Written  in  MatLab® 


%% 

%%  EXTRACT.m 
%% 

%%  Donna  Keating, 

%%  version  24  Feb  96 
%%  29  feb  96 

function  [G]=extract(m) 

global  hdl 
global  img  hdl 
global  imgax 
global  sym 
global  callorient 

if  (nargin  ==  0), 

o/„ - START  UP - 

lastone  =  menu('Most  recent  image  or  new  one??', 'Most  Recent','NEW'); 

if  (lastone  ==  1), 
load  lastdata; 
else 

[filename, pathname]=uigetfile('*  .tif, 'IMAGE  T(X)L',40,40); 
if  filename  ~=  0 

[img,map]=tiffread([pathname, filename]); 
else 

error(['Error  Reading '  filename]); 
end 

end 

save  lastdata 
[rows,cols]=size(img); 

hdl  =  figure('Name',filename,'NumberTitle','ofi',... 

'Units', 'Pixels',... 

'Position', [75  10  max(600,200+cols)  max(600,rows+20)],... 
'Interruptible','yes'); 


imgax  =  gca; 

imghdl  =  imshow(img,map); 

if  (isind(img)),  img  =  ind2gray(img,map);  end 

set(hdl,'UserData'.filename); 
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rtside  =  20+cols+10; 
n  =  fix(cols/5); 

xlims  =  [n  3*n  5*n]; 

set(gca,  'XLimMode', 'manual',... 

'XLim',[l  cols],... 

'YLim',[l  rows],... 

'XTick',xlims,... 

'Units','pixels',... 

'Position', [20  20  cols*l  rows*l],... 
'Intemiptible','yes'); 

cropbutt  =  uicontrol(hdl,'St>'le','Pushbutton',... 

'String', 'SELECT  REGION',... 

'Position', [rtside  140  200  25],... 

'Intemiptible',‘yes',... 

'CallBack',... 

'IMCROPIl;  extract(l) '); 

zoombutt  =  uicontrol(hdl, 'Style', 'Pushbutton',... 
'String','ZOOM',... 

'Position', [rtside  170  200  25],... 

'lnterruptible','yes', . . . 

'CallBack',... 

TMZOOM;  extract(l) '); 

file  =  uimenu(hdl,'Laber,'FILE'); 

menul  =  uimenu(hdl,'Laber,'VIEW'); 

menu4  =  uimenu(hdl,'Label','METRICS'); 

% 

% - 

sym(l)  =  uicontrol(hdl,'Style','radio', ... 

'Position', [rtside  275  150  25], ... 
'String','Horizontal  Bars', ... 

'Value',  1); 

sym(2)  =  uicontrol(hdl,'Style','radio', ... 

'Position', [rtside  300  150  25], ... 
'String','Vertical  Bars', ... 

'Value',0); 

for  i=l:2,  set(sym(i),'UserData',sym(:,[l:(i-l),(i+l);2]));  end 

callorient  =  ['me  =  get(gcf,"CurrentObject");',... 
'if(get(me,"Value")  ==  1),',.,. 
'set(get(me,"UserData"),"Value",0),',... 

'clSG  ' 

'set(me,"Value",l 

'end']; 
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set(sym.'CallBack',callorient); 
% - 


endbutton  =  xucontrol(hdl, 'Style', 'Pushbutton',... 

'String','QUIT',... 

'Position', [rtside  95  150  25],... 
'CallBack','close  all;  clear  all;'); 

newfile  =  uimenu(file,'Laber,'New',... 

'CallBack','cla;  close;  extract'); 

o/„ - .end  start  up - 


elseif  (nargin  ==  1), 


if  (m  ==  1), 


end 

end 

acsl 

return 


function  [bout,rect]  =  imcropii 

%% 

%  IMCROPII  Crop  image. 

%%  IMCROP,  Changed  by  Donna  Keating 
%%  Version  23  Feb  96 

%% 

%%  26  Feb  96,  0710  hrs 

%% 

%%  06  Mar  96,  0500  hrs 

global  hdl 
global  img  hdl 
global  imgax 
global  sidehdl 

%  Determine  which  IMAGE  coordinate  system  is  being  used, 
s  =  [version  '  '];  k  =  find(s<46  &  s>58); 
if  ~isempty(k),  s  =  s(I:min(k));  end 
[ver,count,msg,next]  =  sscanf(s,'%f,l); 

ifver  >  4.1, 
useNew  =  I; 
elseif  ver  <  4.1, 
useNew  =  0; 
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else 

if  s(next)>='a',  useNew  =  1;  else  useNew  =  0;  end 
end 


axes(imgax); 

[x,y,a,hasimage]  =  getimage(imgax); 
if  -hasimage, 

error('The  current  figure  must  contain  an  image  to  use  IMCROP.'); 
end 

rect  =  getrect(hdl)  %  Get  rect  info  from  the  user. 
set(gca,'UserData',rect); 

[m,n]  =  size(a); 

xmin  =  min(min(x));  ymin  =  min(min(y)); 
xmax  =  max(max(x));  ymax  =  max(max(y)); 
if  useNew, 

%  Transform  rectangle  into  row  and  column  indices, 
kx  =  n-1;  ky  =  m-1; 

11  =  max(ceil((rect(2)-ymin)/(ymax-ymin)*ky+l),l); 
jl  =  max(ceil((rect(l)-xmin)/(xmax-xmin)*kx+l),l); 

12  =  min(floor((rect(2)+rect(4)-ymin)/(ymax-ymin)*ky+l),m); 
j2  =  min(floor((rect(l)+rect(3)-xmin)/(xmax-xmin)*kx+l),n); 

else 

%  Transform  rectangle  into  row  and  coliunn  indices.  Fix-up  coordinates 

%  to  deal  with  coordinates  extending  from  pixel  boundaries  rather 

%  than  pixel  centers. 

dx  =  max(  (xmax-xmin)/size(a,2),  eps ); 

dy  =  max(  (ymax-ymin)/size(a,l),  eps ); 

kx  =  (n-l+dx);  ky  =  (m-l+dy); 

11  =  max(ceil((rect(2)-ymin)/(ymax-ymin)*ky+(l-dy/2)),l); 
jl  =  max(ceil((rect(l)-xmin)/(xmax-xmin)*kx+(l-dx/2)),l); 

12  =  min(floor((rect(2)-i-rect(4)-ymin)/(ymax-ymin)*ky+(l-d>72)),ni); 
j2  =  min(fioor((rect(l)+rect(3)-xmin)/(xmax-xmin)*kx+(l-dx/2)),n); 

end 


bout  =  a(il:i2Jl:j2); 


save  img  bout 
sidehdl  =  figure; 

set(sidehdl, 'Position', [3 10  200  400  400]); 

imshow(bout); 

colormap(gray(256)) 

filename  =  get(hdl,'UserData'); 

title(filename) 

clear  bout 

task  =  menu('Bar  Group  OR  Brightness  Panels??','Bar  Group', 'Bright.  Panels'); 
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if  (task  ==  1),  m  =  'fit';  end 
if  (task  ==  2),  m  =  'bp' ;  end 

if(strcmp(m,'fit')), 

goon  =  menu('Analyze  this  bar  group?', 'YES', 'NO'); 
if  (goon  ==  2), 
return 
else 

st=fit(l,0); 

return 

end 

elseif  (strcmp(m,'bp')), 

goon=menu(' Analyze  this  pair  of  Brightness  Panels?', 'YES','NO'); 
if  (goon  ==  2), 
return 
else 

st=fit(l,l); 

return 

end 

else 

return; 

end 


%fit.m 

%% 

%  Version  26  Feb  96,  0151 
%%  12  mar  96  1615 

%%  20  mar  96  0230 

function  [modulation]  =  fit(i,m) 

global  hdl 
global  img  hdl 
global  imgax 
global  sidehdl 
global  sym 
global  callorient 

load  img 

outfile  =  mat2str(['Passltop','.out']); 
figure(sidehdl); 

filename  =  get(hdl,'UserData'); 
colormap(gray(256)) 

img2  =  zeros(size(bout)); 
img2  =  medfilt2(bout); 
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clear  bout 


map  =  gray(256); 
if  (isgray(img2)), 

[imgf,map]  =  gray2ind(img2,256); 
else 

imgf  =  img2; 
end 

clear  img2 

%  Watch  out  for  side  effects  here— am  I  really  changing  img,  or  is  it  pass  by  value? 
if  (m  ==  0), 

%======================================================== 

imshow(imgf,gray(256)); 

title([filename '  Median  Filtered']) 

[M  N]  =  size(img0; 
mask  =  zeros(|M  N]); 
len  =  M*N; 


%■ 


thrl  =  mean(prctile(imgf  ,60)); 
thr2  =  mean(prctile(imgf,60)); 
thr  =  max(thrl,thr2); 

%%  This  is  the  level  I  assume  the  "bar  edge"  is  located  (I  hope  its  valid)  (???) 

imgbw  =  im2bw(imgf,thr); 
imgbw  =  bwmorph(imgbw, 'majority', 2); 
imgbw  =  bwmorph(imgbw,'hbreak',2); 
subplot(  1,2,1) 

imgbw  =  bwperim(imgbw,8); 

imshow(imgbw,2) 

title('Mask') 

%  Use  this  mask  to  find  the  bar  in  the  image,  and 
%  as  a  first  pass  to  approach  a  good  estimate  of  the  modulation 

%  Need  this,  'cause  imgbw  could  be  a  filled  in  square! 

%  Use  imgbw  to  determine  xrange  and  yrange 
hold  on 

%  Want  to  exelude  rows  or  columns  with  fewer  nnz  pixels  than  "some  limiting  value"  that  I 
%  currently  choose  to  be  the  mean  col  or  row  population.  Think  median  will  be  too  high.,.(??) 

xrange  =  fmd(sum(imgbw)  >=  2); 
yrange  =  find(sum(imgbw')>=  2); 

clear  imgbw 


xmax  =  max(xrange); 
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xmin  =  min(xrange); 
ymax  =  max(yrange); 
ymin  =  min(yrange); 
xrange  =  [xmin:xmax]; 
yrange  =  [ymin:ymax]; 

mask(yrange, xrange)  =  ones([length(yrange),length(xrange)]); 
imgff  =  imgf.*mask; 

%% - end  mask - 

XX  =  prctile(imgff(yrange, xrange), 60); 
xy  =  prctile(img£f(yrange, xrange)', 60); 

XX  =  expsmth(xx);  xx  =  expsmth(xx); 
xy  =  expsmth(xy);  xy  =  expsmth(xy); 

%dx2  =  del2(prctile(imgfif(yrange,xrange),60)); 
dx2  =  del2(xx); 

%dy2  =  del2(prctile(imgfiF(yrange,xrange)',60)); 
dy2  =  del2(xy); 

%vardx2  =  var(del2(prctile(inigff(yrange,xrange),60))); 

%vardy2  =  var(del2(prctile(imgff(yrange,xrange)',60))); 
vardx2  =  var(del2(xx)); 
vardy2  =  var(del2(xy)); 

%stddx2  =  std(del2(prctile(img£f(yrange,xrange),60))); 

%stddy2  =  std(del2(prctile(imgff(yrange,xrange)',60))); 
stddx2  =  std(del2(xx)); 
stddy2  =  std(del2(xy)); 

%%  The  DEL2  function  takes  the  "Laplacian"  and  basically  zeros  out 
%%  constant  or  linear  changes!  Making  it  ideal  to  determine  orientation 
%%  With  out  needing  to  know  the  exact  location  of  the  region  of  interest! 

%  Assume  the  orientation  with  largest  variance  is  the  width-wise  direction. 

%  — >If  this  isn't  true,  I  could  fit  a  polynomial  to  the  known  bar  cross-section 
%  and  measme  the  error  of  what  I  have  to  what  I  expect...  noise  could  make 
%  this  unreliable  too-especially  normalizing  the  amplitudes  and  such. 

clear  imgff; 

maskarea  =  bwarea(mask); 
clear  mask; 

meanimg  =  mean(mean(imgf(yrange, xrange))); 

if  (vardx2  >  vardy2), 
orient  =  1; 

grpwidth  =  maskarea/lengthfyrange); 

meanmat  =  (ones([length(yrange),length(xrange)])')*meanimg; 
imgf  =  ((detrend(imgf(yrange, xrange)'))  +  meanmat )'; 
set(sym, 'Value', 0); 
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set(sym,'CallBack',caIlorient); 
elseif  (vardx2  <  vardy2), 
orient  =  -1; 

grpwidth  =  maskarea/length(xrange); 
meanmat  =  ones([length(yrange),length(xrange)])*meanimg; 
imgf  =  detrend(imgf(yrange,xrange))  +  meanmat; 
set(sym, 'Value',  1); 
set(sym,'CalIBack',callorient); 
else 

error('Umesolvable'); 

break; 

end 

fig  =  gcf; 

%  Select  the  best 
if  (orient  ==1), 

crosspit  =  sum(imgf)/length(xrange); 
xvalue=  [l:length(xrange)]; 

betas  =  polyfit(xvalue, crosspit, min(10,length(crossplt)-2)); 
fitplt  =  polyval(betas,xvalue); 
xcen  =  round(sum(crossplt.*xvalue)/sum(crossplt)); 
xcen  =  [xcen-l:xcen+l]; 

cen  =  xcen(find(fitplt(xcen)==max(fitplt(xcen)))) 

cut  =  mean(imgf(:,xcen)');  cut=medfiltl(cut,4); 

df2  =  diff(cut,2); 

sortcut  =  unique(cut); 

upper75  =  round(prctile(sortcut,75)); 

keepers  =  find(cut  >=  upper75); 

lower  =  min(keepers) ;  topper  =  max(keepers); 

yrange  =  lower:topper; 

crosspit  =  sum(imgf)/length(xrange); 

tmpfig  =  figure; 

plot(xvalue,crossplt,'b'); 

[maxloc,minloc,dmaxloc,dminloc,fitplt,mult]  =  finder(crossplt,xvalue,cen); 
hold  on; 

xtra  =  [l:length(fitplt)]/mult; 

plot(xcen,crossplt(xcen),'r*') 

plot(cen,crossplt(cen),'yO') 

plot(xtra,fitplt,'w:'); 

figure(fig); 

cutarea  =  length(yrange)*length(xrange); 
orientlabel  =  'Across  LOF'; 
else 

crosspit  =  sum(imgf  )/length(yrange); 
xx^alue  =  1  :length(yrange); 

betas  =  polyfit(xvalue,crossplt,min(l  l,length(crossplt)-2)); 
fitplt  =  polyval(betas,xvalue); 
ycen  =  round(sum(crossplt.*xvalue)/sum(crossplt)) 
ycen  =  [ycen-l:ycen+l]; 

cen  =  ycen(  find(fitplt(ycen)==max(fitplt(ycen)))) 
cut  =  mean(imgf(ycen,:));  cut  =  medfiltl(cut,4); 
df2  =  diff(cut,2); 
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sortcut  =  imique(cut); 

upper75  =  round(prctile(sortcut,75)); 

keepers  =  find(cut  >=  upper75); 

lower  =  min(keepers) ;  topper  =  max(keepers); 

xrange  =  lower:topper; 

crosspit  =  sum(imgf)/length(yrange); 

tmpfig  =  figure; 

plot(xvalue, crosspit, 'g'); 

[maxloc,minloc,dmaxloc,dininloc,fitplt,mult]=fmder(crossplt,xvalue,cen); 
hold  on; 

xtra  =  [l:length(fitplt)]/mult; 
plot(ycen,crossplt(ycen),'r*') 
plotCcen.crosspltCcen),^©') 
plot(xtra,fitplt,'w: ') ; 
figure(fig); 

cutarea  =  length(yrange)*length(xrange); 
orientlabel  =  'In  LOF'; 
end 

xlabel(orientlabel) 
hold  off; 

fplt  =  fitplt; 
xplt  =  crosspit; 

subplot(l,2,2) 
load  wdgresul 

crosspit  =  ([crossplf  ones([length(crossplt),l])]*LOG_E_GL'); 
fitplt  =  ([fitplt’  ones([length(fitplt),l])]*LOG_E_GL'); 

plotCx^alue,  1 0 ,  ''^crossplt(xvalue),'y') 
hold  on; 

plot(xtra,  1 0 ,  ^fitplt,'w: ') ; 

plot(dmaxloc,10.^crossplt(dmaxloc),'r*'); 

plot(dininloc,10.'''crossplt(dniinloc),'b*'); 

%altitude  =1666;  %  Pass  1,  Average  altitude  in  meters 

%altitude  =  1630;  %  Pass  2 

%altitude  =  1610;  %  Pass  3 

altitude  =  1608;  %  Pass  4 

%altitude  =  1613;  %  Pass  5 

%altitude=  11648; 

%% Which  bar  group?? 

barwidthO  =  grpwidth/6; 
maxloc  =  sort(maxloc); 
minloc  =  sort(minloc); 

if  (length(maxloc)==3  )&(length(minloc)==2), 
barwidthl  =  abs(minloc(l)-maxloc(l)); 
barwidth2  =  abs(maxloc(2)-minloc(l)); 
barwidth3  =  abs(minloc(2)-maxloc(2)); 
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barwidth4  =  abs(maxloc(3)-minloc(2)); 

barwid  =  mean((l/mult)*|barwidthl  barwidth2  barwidth3  barwidth4]) 
barwidth  =  round(barwid) 

stdwid  =  std((l/mult)*Ibarwidthl  barwidth2  barwidth3  barwidth4]); 
else 

goon  =  menu('ERROR;  Not  Resolved', 'Continue??','RETURN'); 
if  (goon  ==  2),  retam;  end 
barwidth  =  barwidthO; 
barwid  =  barwidthO; 
end 

barlength  =  cutarea/grpwidth; 
barspace  =  zeros([l, barwidth]); 
barhat  =  ones([l, barwidth]); 

avgmaxht  =  (mean(crossplt(dmaxloc))); 
avgminht  =  (mean(crossplt(dminloc))); 

deltae  =  abs(avgmaxht-avgminht); 

C  =  lO^deltae; 

K2  =  (C-1)/(C+1) 

avgmaxexp  =  max(10'''abs(avgmaxht),10^abs(avgminht)); 
avgnvinexp  =  niin(10'''abs(avgmaxht),10^abs(avgnunht)); 
meanmod  =  abs(avgmaxexp-avgminexp)/(avgmaxexp+avgminexp); 

[group, con, fwid]  =  fmdgrp(barwid,altitude) 

if  (con  <  0.7), 

guess  =  menu(['Is  this  really  GROUP  #'  num2str(group) '  ?'],'YES','NO'); 
if  (guess  ==  2), 

fixwid  =  menuCCHOOSE  CORRECT  WIDTH',num2str(group+3),nuni2str(group+2),... 

nuin2str(groufrM),niun2str(group-l),num2str(group-2),num2str(group)-3)); 
if  (fixwid  ==  l),group=group+3; 
elseif  (fixwid  ==  l),group=groufr(-3; 
elseif  (fixwid  ==  2),group=groupt-2; 
elseif  (fixwid  ==  3),group=groufrM; 
elseif  (fixwid  ==  4),group=group-l; 
elseif  (fixwid  ==  5),group=group-2; 
elseif  (fixwid  ==  6),group=group-3; 
end 

fwid  =  ( 120E-2/(2^(  l/6))^(3+group))/(altitude*5.737E-6/(3  *2. 54E-2)); 
end 
end 

title(['Group  # '  nuni2str(group)]) 

xlabel(['Width='  nuni2str(fwid)]) 

ylabel('Exposure'); 

bar  =  ones([l, barwidth]); 

barmodel  =ones([  1  ,length(xvalue)])*avgminht; 

model  =  [bar*avgminht  bar*avgmaxht  bar*avgniinht  bar*avgmaxht  bar*avgniinht  bar*avgmaxht 
bar*avgminht]; 


108 


firstmid  =  mult*(  banvid  +  barwid/2  ) 
firstpk  =  maxloc(l) 

err=round(firstpk  -  firstmid); 
if(err>=  1), 

%  model=[ones([l  ,round(err/mult)])*avgminht  model(max(  1  ,abs(err))  :length(model))] ; 

model=model(abs(round(err/mult))+ 1 :  length(model)); 
else 

model  =  model(abs(roimd(err/mult))+l:length(model)); 
end 

barmodel([l;length(model)])  =  model; 
plot(xvalue,10/'barmodel(xvalue),V) 

hold  off 
save  datafit 

filename  =  get(hdl,'UserData'); 
fid  =  fopen(outfile,'a'); 
lprintf(fid,'%s\t  %s\t\n', filename, date); 

params  =  ['group 'fwid 'barwid',  'orient ',  'barwidth ',  'contrast ' 

'K2  ',  'stddel2x ',  'stddel2y  ',  'vardel2x  ' ,  'vardel2y ',  'confidence'  ]; 

fprintf(fid,'%s\n', params); 

format  1  =  ['%i\t  %d\t  %i\t  %i\t  %d\t  %d\t']; 
format2  =  ['%d\t  %g\t  %g\t  %g\t  %g\t  %g\n']; 
params  1  =  [group;fwid;orient;barwidth;barwid;C]; 
params2  =  [K2;stddx2;stddy2;vardx2;vardy2;con]; 
fprintf(fid,format  1 , params  1 ); 
lprintf(fid,format2,params2); 

figure(tmpfig) 

title([filename ',  Group  #',  num2str(group),' , ',  orientlabel]) 

ylabel('Gray  Level') 

xlabel('Pixels') 

drawnow; 

elseif  (m  ==1), 

%========================================================:===== 

imshow(imgf,map); 
pixnum  =  prod(size(imgO); 
title([filename '  Median  Filtered']) 
oldfig  =  gcf; 

xsection  =  medfiltl(mean(imgf),4); 
bins  =  round(2*log(prod(size(imgf)))); 

[xcounts,  xlocate]  =  HlST(xsection,rem(bins,2)+bins); 
sortcounts  =  sort(xcounts); 

L  =  length(sortcounts); 
loci  =  max(find(xcounts  ==  sortcounts(L))); 
loc2  =  max(find(xcounts  ==  sortcounts(L-l))); 
loc3  =  max(find(xcoimts  ==  sortcoimts(L-2))); 
loc4  =  max(find(xcounts  ==  sortcoimts(L-3))); 
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max  I  loc  =xlocate(max(loc  1-1,1):  inin(loc  1+1  ,L)); 
max21oc  =xlocate(max(  1  ,loc2- 1 ) :  min(loc2+ 1  ,L)); 
if  (any(maxlloc==loc3)|any(maxlloc==loc4)|any(maxlloc==loc2)), 
max  1  loc  =  roimd(sum(xcounts(max  1  loc).  *max  1  loc)/sum(xcoimts(max  Hoc))) ; 
else 

maxlloc  =  loci; 
end 

if  (any(max21oc==loc3)|any(max21oc==loc4)|any(max21oc==loc  1 )), 
max21oc  =  round(sxun(xcounts(max21oc).*max21oc)/sum(xcounts(max21oc))); 
else 

max21oc  =  loc2; 
end 

save  xhist  xcounts  xlocate  xsection 
load  wdgresul 

maxloc  =  max(maxlloc,max21oc) 
minloc  =  min(maxlloc,max21oc) 

%%  Assume  histogram  is  counted  at  the  right  edge  of  each  bin,  so  the  mid-value 
%%  must  be  calculated  as  the  average  of  the  selected  bin  and  its  left  neightbor. 
avgmaxht  =  mean(xlocate(maxloc-l:maxloc)); 
avgminht  =  mean(xlocate(minloc-l:  minloc)); 

avgmaxht  =  (abs([avgmaxht  1]*LC)G_E_GL')); 
avgminht  =  (abs((avgminht  l]*LOG_E_GL')); 

delta_e  =  abs(avgmaxht-avgminht); 

C  =  10''delta_e; 

K2  =  (C-1)/(C+1) 
maxht  =  lO^avgmaxht; 
minht  =  lO-^avgminht; 

meanmod  =  abs(maxht-minht)/(maxht+minht); 

save  imghist  L  maxlloc  max21oc  avgmaxht  avgminht 
figure(oldfig); 

filename  =  get(hdl,'UserData'); 
fid  =  fopen(outfile,'a'); 
lprintf(fid,'%s\t  %s\t',filename,date); 

params  =  ['BP '  K2  '  meanmod ', '  contrast ' , '  maxlloc ', '  max21oc ']; 

fprintf(fid,'%s\n',params); 

format  1  =  ['  %d\t  %d\t  %d\t  %d\t  %d\t\n']; 

paramsl  =  [K2;  meanmod;  C;  maxlloc;  max21oc]; 

fprintftfid, format  1  ,params  1 ) ; 

close 

end 

return 


function  [maxloc, minloc, xloc,nloc,fitplt,  L]  =  finder(dataplt,xval,cen) 
%  Finds  Peaks  and  Valleys 


no 


% 

%  Programmer:  Donna  Keating 
% 

%  Search  routine  written  by  Maj  Michael  Roggemann 
%  Modified  by  Doima  Keating 

betas  =  polyfit(xval,dataplt,min(10,length(dataplt)-2)); 
f  =  polyval(betas,xval); 

L=  10; 

x=  [l:length(f)*L]/L; 

cen  =  cen*L; 

fitplt  =  polyval(betas,x); 

fitplt(l:L)=ones([l,L])*fitplt(L+l);  %  Removes  extrapolated  data  less  than  ONE  (causes  false  min). 

In  =  length(fitplt); 
x=  [l:ln]; 

%  Use  pretile  for  initial  values  to  avoid  false  edge  maxima 
%firstmax  =  prctile(fitplt((cen):(ln-l)),30); 
firstmax  =  min(fitplt); 
firstmin  =  prctile(fitplt,60); 

%secondmax  =  prctile(fitplt((cen):(ln-l)),30); 

secondmax  =min(fitplt); 

firstmaxflag  =  0; 

firstminflag  =  0; 

firstmaxloc  =  cen; 

%find  min  and  max  going  right 

for  j=(cen-l):(ln-l), 

%  CASE  1 :  The  next  element  is  larger,  but  max  not  yet  found 

if  (fitpltfj)  >=  firstmax)  &  (firstmaxflag  ==  0), 

firstmax  =  fitplt(j); 
firstmaxloc  =  j; 

if  (fitplt(j+l)  <=  fitpltCj)), 

firstmaxflag  =  1; 

end; 

end; 

%  Found  first  max,  now  need  to  find  min! 
if  (flrstmaxflag==l)  &  (firstminflag==0), 

if  (fitpltO)  <=  firstmin  ), 
firstmin  =  fitplt(j); 
firstminloc  =  j; 
end; 
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if(fitplt(j+l)>=fitplt(j))&(fitplt(j-2)>=fitplt(j)), 
firstminflag  =  1; 

end; 

end; 

if  (firstmaxflag==l)  &  (firstminflag==l)  &  (fitplt(j)>=secondmax), 

if  (fitplt(j)>=fitplt(j-l))&(fitplt(j)>=fitplt(j+2)), 
secondmax  =  fitplt(j); 
secondmaxloc  =  j; 

elseif  (j==(ln-l))  &  (fitplt(j+l)>fitplt(j)), 
secondmax  =  fitplt(j+l); 
secondmaxloc  =  j+l  ; 
end; 
end; 

end; 

%fmd  min  and  max  going  left 

thirdmax  =  prctile(fitplt,30); 
secondmin  =  min(firstmax,  secondmax); 
thirdmaxflag  =  0; 
secondminflag  =  0; 

%fmd  min  and  max  going  right 

for  j=firstmaxloc :  - 1 :  (L+ 1 ), 

if  (thirdmaxflag  ==  0)  &  (secondminflag  ==0), 

if  (fitplt(j)<=secondmin), 
secondmin  =  fitplt(j); 
secondminloc  =  j; 
end; 

if  (fitplt(j-l)>=fitplt(j))&(fitplt(j+2)>=fitplt(j)), 
secondminflag  =  1 ; 
end; 

end; 

if  (fitplt(j)  >=  thirdmax)  &  (thirdmaxflag  ==  0)  &  (secondminflag  ==  1), 
thirdmax  =  fitplt(j); 
thirdmaxloc  =  j; 

if  (fitplt(j-l)  <=  fitplt(j))  &  (fitplt(j+l)  <=  fitpltfj)), 
thirdmaxflag  =  1; 
end; 
end; 
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end; 


%-■ 


maxloc=[firstmaxloc  secondmaxloc  thirdmaxloc]; 
minloc=[firstminloc  secondminloc]; 

%  Convert  from  fitplt  coordinates  back  to  dataplt  coords, 
xloc  =  unique(maxloc)/L; 
nloc  =  unique(ininloc)/L; 

for  i=  1  :length(xloc), 
lox  =  floor(xloc(i)); 
hix  =  ceil(xloc(i)); 
xbounds  =  [lox:hix]; 
if  (dataplt(lox)>=dataplt(hix)), 

X  =  lox; 
else 

X  =  hix; 
end 

xloc(i)  =  X; 
end 

fori  =  l:length(nloc), 
lox  =  floor(nloc(i)); 
hix  =  ceil(nloc(i)); 
xbounds  =  [lox:hix]; 
if  (dataplt(lox)<=dataplt(hix)), 

X  =  lox; 
else 

X  =  hix; 
end 

nloc(i)  =  X; 
end 

save  finder 


%%  Findgrp.m 

function  [group, confid,twid]  =  findgrp(barwidth, altitude) 

%%  ASSUME  Sensor  is  the  KS-87 

%%  I  have  pixels  of  bar  width  and  ground  meters  of  bar  width  — »  Which  bar 
%%  Need  ground  width  per  pixel! !  I  Need  ground  pitch  (how  much  ground  distance  in  a 
%%  5.737  um  image  plane  distance??? 

%%  pitch:FL,  as  ground  pitchialtitude 

pitch  =  5.737*10''(-6); 
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%  meters 

FL  =  3*2.54*10^(-2); 

%FL=  18*2.54*10^(-2); 

%  meters 

groundpitch  =  (pitch/FL)*altitude; 

%  So,  estimated  width  of  the  bar  in  question  is  groundpitch*pixelmunber! ! ! 
gbarwidth  =  barwidth*groundpitch 
%  So  in  what  known  bargroup  range  does  this  fall???? 
prog=  2''(l/6); 

%  Bar  dimension  reduction  factor, 
biggest  =  120*10'^(-2); 

%  meters 

if  (gbarwidth  >=  biggest/(prog'^4)), 
left  =  biggest-gbarwidth; 
right  =  gbarwidth-biggest/prog'^4; 
if  (left  <  right), 
group  =  0; 

confid =1  -  left/right; 
twid  =  biggest/groundpitch; 
return; 
else 

group  =  1; 

confid=  1  -  right/(left); 
twid  =  (biggest/(progM))/groundpitch; 
return; 
end 
end 

for  i=  1:12, 

leftbar  =  biggest/prog^(3+i); 
rightbar  =  biggest/prog^(3+i+l); 
if  (gbarwidth<=leftbar  &  gbarwidth>=rightbar) 
left  =  leftbar-gbarwidth; 
right  =  gbarwidth-rightbar; 
delta  =  left-right; 
if  (left  <  right), 
group  =  i; 

confid  =  1  -  left/right; 
twid  =  leftbar/groundpitch; 
return; 
else 

group  =  i+l; 
confid  =  1  -  right/left; 
twid  =  rightbar/groundpitch; 
return; 
end 
end 
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% - 

end  %forloop 


% 

%  expsmth.m,  exponentially  weighted  smoothing 
% 

function  result  =  expsmth(DATA) 

ALPHA  =  0.2; 

f  =  zeros([length(DATA)+l,l]); 
result  =  zeros([length(DATA)+l,l]); 

f(l)  =  minCDATA); 

f(2)  =  DATA(1)*ALPHA  +  (l-ALPHA)*f(l); 
for  i  =  l:length(DATA), 

f(i+2)  =  DATA(i)*ALPHA  +  (l-ALPHA)*f(i+l); 
end; 

result  =  f; 


function  y  =  prctile(x,p); 

%PRCTILE  gives  the  percentiles  of  the  sample  in  X. 

%  Y  =  PRCTILE(X,P)  returns  a  value  that  is  greater  than  P  percent 
%  of  the  values  in  X.  For  example,  if  P  =  50  Y  is  the  median  of  X. 

% 

%  P  may  be  either  a  scalar  or  a  vector.  For  scalar  P,  Y  is  a  row 

%  vector  containing  Pth  percentile  of  each  column  of  X.  For  vector  P, 

%  the  ith  row  of  Y  is  the  P(i)  percentile  of  each  column  of  X. 

%  Copyright  (c)  1993  by  The  MathWorks,  Inc. 

%  $Revision;  1.1  $  $Date:  1993/05/24  18:56:10  $ 


[prows  pcols]  =  size(p); 
if  prows  ~  1  &  pcols  —  1 
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error('P  must  be  a  scalar  or  a  vector.'); 
end 

if  any(p  >  100)  |  any(p  <  0) 

error('P  must  take  values  between  0  and  100'); 
end 

XX  =  sort(x); 

[m,n]  =  size(x); 

ifm==l  |n==l 
m  =  max(m,n); 
n=  1; 

q  =  100*(0.5:m  -  0.5)./m; 

XX  =  [min(x);  xx(:);  max(x)]; 
else 

q=  100*(0.5:m  -  0.5)./m; 

XX  =  [min(x);  xx;  max(x)]; 
end 


q=[0ql00]; 


y  =  interpl(q,xx,p); 
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